load packages:

library(googlesheets4) #scrapes google sheets
library(dplyr) # dataframe manipulation
library(ggplot2) # plotting package
library(lubridate)
library(tidyr)
library(readr)
library(scales)
library(growthrates) 
library(patchwork)
library(grid)
library(qiime2R)
library(phyloseq)
library(viridis)
library(magrittr)
library(vegan)
library(stringr)
library(ape)

Create a color pallete

pp_palette <- c(
  # Original + Extended Palette (50 colors)
  "#7D8F69", "cornflowerblue", "#F5C45E", "#C1D3D5", "plum",
  "hotpink4", "lightsalmon", "steelblue", "indianred", "springgreen4", 
   "#102E50", "maroon", "lightgreen", "lightblue", "gold",
  "#BFA5A0", "#A26769", "blue", "red", "lightpink", 
  "#71816D", "#FFE6E6", "#E5C1CD", "#A3B18A", "salmon", "magenta", "orchid","deepskyblue", "tomato",
    "cyan", "hotpink", "yellow", "orange",
   "#C0A98E", "#7A9E9F", "#A1869E", "#BFD8B8", "#F7DD72","green", "darkgreen","purple", 
  "#FFE5AD",
  "#CDB4DB"  
)
scales::show_col(pp_palette)

1 Amplicons associated with P. bahamense isolates from Indian River Lagoon (IRL) (FWC1003) and Old Tampa Bay (OTB) FLorida (FWC1006)

Import sequence data from QIIME 2: Amplicons were processed on a QIIME 2 pipeline, which includes DADA denoising. Seqeunces were classified with Native Bayers Classifier

ps<-qza_to_phyloseq(features="table.qza")

Import and prepare metadata:

meta <- read.csv("SampleMetaData.csv")
meta$SampleID <- paste0(meta$SampleID, "_6283")
row.names(meta)<- meta$SampleID
meta$TreatRep <- paste(meta$PCRDate_Treat, meta$SampleInfo_Rep,  sep = "_")
#keep only the P. bahamense project
pyro_only <- meta[grepl("^Pyro", meta$Project, ignore.case = TRUE), ]
META <- sample_data(pyro_only)

Import and prepare taxonomy data:

taxonomy <- read.csv("taxonomy.csv", stringsAsFactors = FALSE)

names(taxonomy) <- c("row", "tax", "Confidence") 
row.names(taxonomy) <-taxonomy[[1]] #move the feature ID column to become the row names
taxonomy <- taxonomy[,(-1)] #delete the feature ID  column 
taxonomy <-  separate(taxonomy, tax, c("Domain","Phylum", "Class", "Order", "Family", "Genus", "Species", "D7", "D8", "D9", "D10", "D11", "D12", "D13", "D14"), sep = ";", fill = "right")
taxonomy <- taxonomy[,c(1:7)] #keep first seven levels

taxonomy$D0 <- with(taxonomy, ifelse(Order == " o__Chloroplast", "Chloroplast", "Bacteria")) #create column for ordering

# Reorder columns and turn into phyloseq tax_table
col_order<- c("D0", "Domain","Phylum", "Class", "Order", "Family", "Genus", "Species")
taxonomy<- taxonomy[, col_order]
taxmat <- as.matrix(taxonomy)
TAX = tax_table(taxmat)

Subsetting the phyloseq object:

#Keep only Pyro samples
ps <- prune_samples(rownames(pyro_only), ps)

#Remove OTUs that have zero counts in all Pyro samples
ps <- prune_taxa(taxa_sums(ps) > 0, ps)
ps <- merge_phyloseq(ps, TAX, META)

Bar plot at domain level:

#Keep only Bacteria + Chloroplast
ps<- subset_taxa(ps, D0 %in% c("Bacteria", "Chloroplast") )

#Convert counts to relative abundance (%)
psra<- transform_sample_counts(ps, function(OTU) 100* OTU/sum(OTU)) 

#Group OTU by DO level
glomD1<- tax_glom(psra, 'D0')

#plot RA
taxabarplot <- plot_bar(glomD1, x = "SampleID", fill = "D0") +
  scale_y_continuous(expand = c(0, 0)) +
  theme_classic() +
  theme(
    legend.title = element_blank(),
    text = element_text(size = 14),
    axis.text.x = element_text(angle = 90)
  ) +
  ylab("Relative Abundance (%)") +
  xlab("") +
  scale_fill_manual(values = c("lightgrey", "#7BB03B")) +
  facet_grid(~Project, scales = "free_x")


taxabarplot

Remove chloroplast sequences from data:

ps_nochloro <- subset_taxa(ps, Order != " o__Chloroplast" & Domain != "Unassigned")

Calculate observed richness:

plugin <- ps_nochloro %>%
            estimate_richness(measures = "Observed") %$% Observed

#Pull project metadata field from phyloseq object
Project <- ps_nochloro %>% sample_data %$% Project

#Combine richness + project into one dataframe
richness<- data.frame(plugin, Project)
names(richness) <- c("richness", "Project")


richness %>%group_by(Project) %>% summarize(mean = mean(richness), min = min(richness), max = max(richness))
## # A tibble: 1 × 4
##   Project  mean   min   max
##   <chr>   <dbl> <dbl> <dbl>
## 1 Pyro     97.5    57   295

Relative abundance and sums the counts of all OTUs that belong to the same Order:

ps_nochloro_RA<- transform_sample_counts(ps_nochloro, function(OTU) 100* OTU/sum(OTU))

ps_nochloro_RA_glomO<- tax_glom(ps_nochloro_RA, 'Order')

Bar plot at the order level:

# subset by pyro samples
pyro <- subset_samples(ps_nochloro_RA_glomO, Project == "Pyro")

# Remove low abundance orders
pyro = filter_taxa(pyro, function(x) sum(x) > 1, TRUE)
                       
taxabarplot <- plot_bar(pyro, x= "SampleID", fill = "Order") +  scale_y_continuous(expand = c(0, 0)) + ggtitle("")  + theme(legend.title=element_blank()) + geom_bar(aes(fill=Order), stat="identity", position="stack", width =0.9) +theme_classic() + theme(text = element_text(size=14))+theme(axis.text.x = element_text(angle = 90)) + xlab("Sample") + ylab("Relative Abundance (%)") + theme(text = element_text(size=14)) + scale_fill_manual(values= pp_palette)+ xlab("")

taxabarplot

Bar plot at family level:

ps_nochloro_RA<- transform_sample_counts(ps_nochloro, function(OTU) 100* OTU/sum(OTU))

#Group OTU by family and sum abundances
ps_nochloro_RA_glomF<- tax_glom(ps_nochloro_RA, 'Family')

#Keep only the Pyro samples
pyro <- subset_samples(ps_nochloro_RA_glomF, Project == "Pyro")

#Remove low abundance orders
pyro = filter_taxa(pyro, function(x) sum(x) > 1, TRUE)
                       
taxabarplot<-plot_bar(pyro, x= "SampleID", fill = "Family") +  scale_y_continuous(expand = c(0, 0)) + ggtitle("")  + theme(legend.title=element_blank()) + geom_bar(aes(fill=Family), stat="identity", position="stack", width =0.9) +theme_classic() + theme(text = element_text(size=14))+theme(axis.text.x = element_text(angle = 90)) + xlab("Sample") + ylab("Relative Abundance (%)") + theme(text = element_text(size=14)) + scale_fill_manual(values= pp_palette)+ xlab("")

taxabarplot

Parse data:

#Filter SampleID
pyroE <- subset_samples(pyro, !(SampleID %in% c("L14_6283", "L15_6283", "L17_6283", "L18_6283")))

#Remove low abundance orders
pyroE = filter_taxa(pyroE, function(x) sum(x) > 1, TRUE)
                       
taxabarplot<-plot_bar(pyroE, x= "SampleID", fill = "Family") +  scale_y_continuous(expand = c(0, 0)) + ggtitle("")  + theme(legend.title=element_blank()) + geom_bar(aes(fill=Family), stat="identity", position="stack", width =0.9) +theme_classic() + theme(text = element_text(size=14))+theme(axis.text.x = element_text(angle = 90)) + xlab("Sample") + ylab("Relative Abundance (%)") + theme(text = element_text(size=14)) + scale_fill_manual(values= pp_palette)+ xlab("")

taxabarplot

Determine Unique OTU’s and abundances in P. bahamense data:

pE_m <- psmelt(pyroE)
unique(pE_m$OTU)
## [1] "12a986d7ed4f3be699daa0feaba41d8a" "07ae9a3ebf44d4a821d001debe7d9b0c"
## [3] "e9034a72595933ebd5c5df08eb4a08fc"
  1. 12a986d7ed4f3be699daa0feaba41d8a o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

12a986d7ed4f3be699daa0feaba41d8a GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGCGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

uncultured sanger sequence from Mangrove Sediment in Hong Kong - GenBank: MH091208.1, MH091199.1 (2 best hits, 100% ) Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs

Diversity of salt marsh prokaryotes, M.A. Moran, uncultured Hyphomicrobiaceae bacterium - lon=81.2797W, lat=31.3884N; sediment 14-16cm collected on Feb 01, 2002, Sapelo Island Microbial Observatory Dean Creek Marsh sampling site” (3rd best hit 99%)

pE_m %>%  filter(OTU == "12a986d7ed4f3be699daa0feaba41d8a") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##       mean      min      max
## 1 40.82156 27.43475 50.63596
  1. e0169d54945a8a584037e96365ec7bb3 o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

e0169d54945a8a584037e96365ec7bb3 GCTGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGCGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

MH091208.1 uncultured bacteria - Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs - isolation_source=“mangrove sediments” in Hong Kong, 99%

pE_m %>%  filter(OTU == "e0169d54945a8a584037e96365ec7bb3") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##   mean min  max
## 1  NaN Inf -Inf
  1. aa0c627da86b1ddbba980c9dfd6ac380 o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

aa0c627da86b1ddbba980c9dfd6ac380 GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGCGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGACTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

MH091208.1 uncultured bacteria - Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs - Hong Kong isolation_source=“mangrove sediments” , 99%

pE_m %>%  filter(OTU == "aa0c627da86b1ddbba980c9dfd6ac380") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##   mean min  max
## 1  NaN Inf -Inf
  1. f6a0f1f2d2a083164d1fb145b000b96b

o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

f6a0f1f2d2a083164d1fb145b000b96b GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGTGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

MH091208.1 Uncultured bacteria Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs - isolation_source=“mangrove sediments” - 99%

pE_m %>%  filter(OTU == "f6a0f1f2d2a083164d1fb145b000b96b") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##   mean min  max
## 1  NaN Inf -Inf
  1. ff08ec2d8a67e142bc459567275bb888 o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

ff08ec2d8a67e142bc459567275bb888 GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGCGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGACTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

MH091208.1 uncultured bacteria Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs Hong Kong isolation_source=“mangrove sediments” 99%

EU488009.1 uncultured bacteria Characterization of the lucinid bivalve-bacteria symbiotic system: the significance of the geochemical habitat on bacterial symbiont diversity and phylogeny isolation_source=“siliciclastic sedment from Thalassia sea grass bed” 99%

other marine sediments, estuary sediments

pE_m %>%  filter(OTU == "ff08ec2d8a67e142bc459567275bb888") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##   mean min  max
## 1  NaN Inf -Inf

  1. 07ae9a3ebf44d4a821d001debe7d9b0c
    c__Acidimicrobiia o__Microtrichales f__uncultured g__uncultured NA

07ae9a3ebf44d4a821d001debe7d9b0c GCAGCAGTGGGGAATCTTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGCGGGAAGACGGCCTTCGGGTTGTAAACCGCTTTCAGCAGGGACGAAATTGACGGTACCTGCAGAAGAAGCTCCGGCCAACTACGTGCCAGCAGCCGCGGTAAGACGTAGGGGGCGAGCGTTGTCCGGAATCATTGGGCGTAAAGGGCTCCTAGGTGGTTCAGTAAGTCGACTGTGAAAATCCAAGGCTCAACCTTGGGACGCCAGTCGATACTGCTGTGACTCGAGTTCGGTAGAGGAGTGTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCAACGGCGAAGGCAGCACTCTGGGCCGATACTGACACTGAAGAGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

KM840989.1 uncultured bacteria - Microbial diversity of indigenous bacteria in a 129I contaminated groundwater plume at the Hanford Site, Washington - Sanger, isolation_source=“iodine (I-129) contaminated groundwater” 100%

HF558551.1 uncultured bacteria - Iron- and Sulphur- cycling bacteria mobilize copper in a multiple extreme mine tailings in the Atacama Desert, Chile - isolation_source=“tailing material” - 100%

JQ427269.1 uncultured bacteria - Bacterial diversity in an alkaline saline soil spiked with anthracene, samger sequenced, isolation_source=“soil”, 100%

JQ665390.1 uncultured actinobacterium, Diversity of unculturable Actinomycetes in coastal wetlands of the Yellow River estuary, isolation_source=“soil sample from coastal wetlands”, 99%

pE_m %>%  filter(OTU == "07ae9a3ebf44d4a821d001debe7d9b0c") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##       mean      min      max
## 1 33.68851 22.47499 44.42002
  1. 9b8accad19c30d7a69a6290553111166 c__Acidimicrobiia o__Microtrichales f__uncultured g__uncultured NA

9b8accad19c30d7a69a6290553111166 GCTGCAGTGGGGAATCTTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGCGGGAAGACGGCCTTCGGGTTGTAAACCGCTTTCAGCAGGGACGAAATTGACGGTACCTGCAGAAGAAGCTCCGGCCAACTACGTGCCAGCAGCCGCGGTAAGACGTAGGGGGCGAGCGTTGTCCGGAATCATTGGGCGTAAAGGGCTCCTAGGTGGTTCAGTAAGTCGACTGTGAAAATCCAAGGCTCAACCTTGGGACGCCAGTCGATACTGCTGTGACTCGAGTTCGGTAGAGGAGTGTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCAACGGCGAAGGCAGCACTCTGGGCCGATACTGACACTGAAGAGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

KM840989.1 uncultured bacteria - Microbial diversity of indigenous bacteria in a 129I contaminated groundwater plume at the Hanford Site, Washington - isolation_source=“iodine (I-129) contaminated groundwater” - 99%

pE_m %>%  filter(OTU == "9b8accad19c30d7a69a6290553111166") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##   mean min  max
## 1  NaN Inf -Inf

  1. e9034a72595933ebd5c5df08eb4a08fc o__Oceanospirillales f__Alcanivoracaceae1 g__Alcanivorax s__Alcanivorax_venustensis

e9034a72595933ebd5c5df08eb4a08fc GCAGCAGTGGGGAATCTTGGACAATGGGGGCAACCCTGATCCAGCCATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTCAGCAGGGAGGAAGGCTTACCCCTAATACGGGTGAGTACTTGACGTTACCTGCAGAAGAAGCACCGGCTAATTTCGTGCCAGCAGCCGCGGTAATACGAAAGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGCGGTGTGTTAAGTCGGATGTGAAAGCCCAGGGCTCAACCTTGGAATTGCATCCGATACTGGCACGCTAGAGTGCAGTAGAGGGAGGTGGAATTTCCGGTGTAGCGGTGAAATGCGTAGAGATCGGAAGGAACACCAGTGGCGAAGGCGGCCTCCTGGACTGACACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

PP516259.1 organism=“Alloalcanivorax venustensis” Exploration and conservation of bacterial community from the Arabian Sea seamount isolation_source=“Arabian Sea water”, 100%

pE_m %>%  filter(OTU == "e9034a72595933ebd5c5df08eb4a08fc") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##      mean      min      max
## 1 25.4058 18.12139 34.43303

  1. c02844d75d512c1510de7334b6687731 o__Rhizobiales f__Hyphomicrobiaceae g__uncultured s__uncultured_bacterium

c02844d75d512c1510de7334b6687731 GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAGCTCTTTTGGCGGGGAAGATAATGACGGTACCCGCAGAATAAGCTCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGAGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGATTTGTTAGTCAGGGGTGAAATCCCGGGGCTCAACCCCGGAACTGCCTTTGATACTGCAAATCTCGAGTCCGAGAGAGGTGGGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCGGTGGCGAAGGCGGCCCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

OK235762.1 Uncultured Pedomicrobium sp. isolation_source=“date palm rhizosphere soil” Saudi Arabia - 99%

HQ697761.1 Uncultured bacteria isolation_source=“hydrocarbon contaminated saline-alkali soil” China - 99%

99% to several other contaminated soil sites (real and experimental)

KP098952.1 Uncultured bacteria - The microbiome of methanol-utilizing denitrification systems contains new bacterial groups - isolation_source=“denitrification bioreactor” - 99%

pE_m %>%  filter(OTU == "c02844d75d512c1510de7334b6687731") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
##   mean min  max
## 1  NaN Inf -Inf

Plot shannon diversity and relative abundance:

# -------------------------------
# Subset Pyro samples and filter
# -------------------------------
pyro <- subset_samples(ps_nochloro_RA_glomO, Project == "Pyro")
pyro <- subset_samples(pyro, Treatment != "AB")
pyro <- filter_taxa(pyro, function(x) sum(x) > 1, TRUE)
pyro <- subset_samples(pyro, 
                       Treatment %in% c("1003", "Replete") & 
                       TreatRep %in% c("1003_1", "1_2"))

# Set factor levels for Notes_Filter
pyro@sam_data$Notes_Filter <- factor(as.character(pyro@sam_data$Notes_Filter),
                                     levels = c("10", "2"))

# Clean Order names in tax_table
tax_table(pyro)[, "Order"] <- gsub(" o__", "", tax_table(pyro)[, "Order"])

# -------------------------------
# Plot relative abundance barplot
# -------------------------------
barplot <- plot_bar(pyro, x = "TreatRep", fill = "Order") +
  geom_bar(aes(fill = Order), stat = "identity", position = "stack", width = 1) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_manual(values = pp_palette) +
  facet_grid(Notes_Filter ~ Treatment, scales = "free") +
  scale_x_discrete(
    expand = c(0.5, 0.05),
    labels = c(
      "1003_1" = "IRL",
      "1_2"    = "OTB"
    )
  ) +
  theme_classic() +
  theme(
    legend.title = element_blank(),
    legend.position = "right",
    legend.justification = "center",
    legend.box = "horizontal",
    axis.text.x = element_text(angle = 0, vjust = 0.5, size = 16, color = "black"),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.title.x = element_text(size = 20, color = "black"),
    axis.title.y = element_text(size = 20, color = "black"),
    axis.ticks.x = element_line(color = "black", linewidth = 0.9),
    axis.ticks.y = element_line(color = "black", linewidth = 0.9),
    strip.text.x = element_blank(),  # remove top facet labels
    strip.text.y = element_text(size = 20, color = "black"), # row facets
    text = element_text(size = 20, color = "black")
  ) +
  ylab(expression("Relative Abundance (%)")) +
  xlab("")

barplot

# -------------------------------
# Calculate Shannon index
# -------------------------------
pyro <- subset_samples(ps_nochloro_RA_glomO, Project == "Pyro")
pyro <- subset_samples(pyro, TreatRep %in% c("1003_1", "1_2", "1_1", "1_3"))
pyro <- subset_samples(pyro, Treatment != "AB")

shannon_df <- estimate_richness(pyro, measures = "Shannon")

# Extract metadata and combine with diversity
sample_meta <- as(sample_data(pyro), "data.frame")
shannon_df <- cbind(shannon_df, sample_meta)

# -------------------------------
# Plot Shannon diversity
# -------------------------------
shannon_plot <- ggplot(shannon_df, aes(x = Treatment, y = Shannon, color = Notes_Filter)) +
  
  # Color mapping
  scale_color_manual(values = c("10" = "navy", "2" = "maroon")) +
  
  # Points
  geom_point(size = 5) +
  
  # Custom x-axis labels
  scale_x_discrete(labels = c(
    "1003" = "IRL",
    "Replete" = "OTB"
  )) +
  
  # Classic theme
  theme_classic() +
  
  # Axis labels
  ylab("Shannon Diversity Index") +
  xlab("") +
  labs(color = "Filter Size") +  

  # Theme adjustments
  theme(
    legend.title = element_text(size = 20),
    axis.text.x = element_text(angle = 0, vjust = 0.5, size = 16, color = "black"),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.title.x = element_text(size = 20, color = "black"),
    axis.title.y = element_text(size = 20, color = "black"),
    axis.ticks.x = element_line(color = "black", linewidth = 0.9),
    axis.ticks.y = element_line(color = "black", linewidth = 0.9),
    strip.text.x = element_blank(),  # remove top facet labels
    strip.text.y = element_text(size = 20, color = "black"), # row facets
    text = element_text(size = 20, color = "black"))

# Print Shannon plot
shannon_plot

1.1 Figure 2

bo <- shannon_plot + barplot
bo

2 Growth of a Low-Diversity P. bahamense Isolate

Chlorophyll-a fluorescence (RFU) data from P. bahamense cultures were used to visualize growth patterns, quantify maximum biomass, and estimate specific growth rates under media conditions deficient in media lacking one targeted B vitamin. Specific growth rates were computed using the “alleasy_linear” function of the “growthrates” package in R Studio (Hall et al., 2014; RStudio Team, 2020).

Import the Dataset:

fluor<-as.data.frame(read_csv("fluor.csv"))
fluor$Time <- as.numeric(sub('.', '', fluor$Time))
fluor$Date <- mdy(fluor$Date)
fluor$Time <- fluor$Time * 48 /24

2.1 Visulaize growth:

LD_OTB_all_growth <- ggplot(
 fluor %>%
    filter(
      Treatment %in% c("-Cobalamin", "Replete", "-Thiamine", "-Biotin"),
      Round %in% c(1, 2,3,4)),
  aes(x = Time, y = `Chlorophyll-A (RFU)`, color = Treatment)) +
  geom_point(size = 3)  + geom_smooth(se=FALSE)+
  # Manual color scale
  scale_color_manual(
    values = c(
      "-Cobalamin" = "#102E50",
      "Replete"    = "#F5C45E",
      "-Thiamine"  = "deeppink4",
      "-Biotin"    = "#667028"
    ),
    labels = c(
      "-Cobalamin" = expression(" —B"[12]),
      "-Thiamine"  = expression("—thiamine"),
      "-Biotin"    = expression("—biotin"),
      "Replete"    = "replete"),
    name = "Treatment") +
 labs(
    x = "Time (days)",
    y = expression(paste("Chlorophyll-", italic("a"), " fluorescence (RFU)"))
  ) +
  facet_grid(~Round) +
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.title = element_text(size = 20),
    axis.text  = element_text(size = 20, color = "black"),
    strip.text = element_text(size = 20, color = "black"),
    legend.title = element_text(size = 18, color = "black"),
    legend.text  = element_text(size = 18, color = "black"),
    legend.position = "bottom",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  )

LD_OTB_all_growth

2.2 Calculate maximum RFU values for each round and medium treatment

max_bio <- fluor %>%
  filter(Round %in% c(1, 2, 3, 4)) %>%
  filter(Treatment %in% c("-Cobalamin", "Replete", "-Biotin", "-Thiamine")) %>%
  group_by(Treatment, Round, Rep) %>%
  summarise(Max_Chl = max(`Chlorophyll-A (RFU)`, na.rm = TRUE), .groups = "drop")

print(max_bio)
## # A tibble: 48 × 4
##    Treatment Round   Rep Max_Chl
##    <chr>     <dbl> <dbl>   <dbl>
##  1 -Biotin       1     1   2911.
##  2 -Biotin       1     2   3220.
##  3 -Biotin       1     3   6239.
##  4 -Biotin       2     1   6771.
##  5 -Biotin       2     2   4230.
##  6 -Biotin       2     3   3492.
##  7 -Biotin       3     1   4288.
##  8 -Biotin       3     2   3237.
##  9 -Biotin       3     3   1767.
## 10 -Biotin       4     1   2049.
## # ℹ 38 more rows

Run an ANOVA on maximum RFU values and Tukey HSD post-hoc test for treatment

aov_max_bio <-  aov(Max_Chl ~ Treatment, data = max_bio)
TukeyHSD(aov_max_bio)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Max_Chl ~ Treatment, data = max_bio)
## 
## $Treatment
##                           diff       lwr       upr     p adj
## -Cobalamin--Biotin   -3365.871 -4602.055 -2129.687 0.0000000
## -Thiamine--Biotin     1961.220   725.036  3197.404 0.0006407
## Replete--Biotin       -757.515 -1993.699   478.669 0.3694455
## -Thiamine--Cobalamin  5327.091  4090.907  6563.275 0.0000000
## Replete--Cobalamin    2608.356  1372.172  3844.540 0.0000068
## Replete--Thiamine    -2718.735 -3954.919 -1482.551 0.0000030
summary(aov_max_bio)
##             Df    Sum Sq  Mean Sq F value   Pr(>F)    
## Treatment    3 174966709 58322236   45.35 1.63e-13 ***
## Residuals   44  56590796  1286154                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prepare data for plotting maximum RFU values:

# Calculate mean and SD
summary_stats <- max_bio %>%
  group_by(Treatment, Round) %>%
  summarise(
    mean_chl = mean(Max_Chl, na.rm = TRUE),
    sd_chl   = sd(Max_Chl, na.rm = TRUE),
    .groups = "drop"
  )

# Make Treatment a factor in both datasets
max_bio$Treatment <- factor(max_bio$Treatment, 
                            levels = c( "-Biotin", "-Cobalamin","-Thiamine", "Replete"))

summary_stats$Treatment <- factor(summary_stats$Treatment, 
                                  levels = c("-Cobalamin", "-Thiamine", "-Biotin", "Replete"))

max_bio$Treatment <- factor(
  max_bio$Treatment,
  levels = c("-Cobalamin", "Replete", "-Thiamine", "-Biotin")
)

summary_stats$Treatment <- factor(
  summary_stats$Treatment,
  levels = c("-Cobalamin", "Replete", "-Thiamine", "-Biotin"))

Plot maximum RFU values:

max_bio_LD <- ggplot(max_bio, aes(x = Treatment, y = Max_Chl, color = Treatment)) +
  
  # Raw points (per round)
  geom_point(size = 3, shape = 16, position = position_jitter(width = 0.1)) +
  
  # Error bars: mean ± SD across replicates
  geom_errorbar(
    data = summary_stats,
    aes(
      x = Treatment,
      y = mean_chl,
      ymin = mean_chl - sd_chl,
      ymax = mean_chl + sd_chl,
      color = Treatment
    ),
    width = 0.3,
    position = position_dodge(width = 0.6),
    inherit.aes = FALSE
  ) +
  
  # Open circles for replicate-averaged means
  geom_point(
    data = summary_stats,
    aes(x = Treatment, y = mean_chl, color = Treatment),
    shape = 21,
    size = 5,
    stroke = 1.2,
    position = position_dodge(width = 0.6),
    inherit.aes = FALSE
  ) +
  
  # Labels
  labs(
    x = "Media Treatment",
    y = expression(paste("Maximum RFU"))
  ) +
  
  # Manual colors
  scale_color_manual(
    values = c(
      "-Cobalamin" = "#102E50",
      "-Thiamine"  = "deeppink4",
      "-Biotin"    = "#667028",
      "Replete"    = "#F5C45E"
    ),
    labels = c(
      "-Cobalamin" = expression("—B"[12]),
      "-Thiamine"  = expression("—Thiamine"),
      "-Biotin"    = expression("—Biotin"),
      "Replete"    = "Replete"
    ),
    name = "Treatment",
    guide = guide_legend(
      override.aes = list(
        shape = 21,
        fill  = "white",
        stroke = 1.2,
        size = 2
      )
    )
  ) +
  
  # X-axis label formatting
  scale_x_discrete(
    labels = c(
      "-Biotin"    = expression("\u2013Biotin"),
      "-Cobalamin" = expression("\u2013B"[12]),
      "-Thiamine"  = expression("\u2013Thiamine"),
      "Replete"    = "Replete"
    )
  ) +
  
  scale_y_continuous(labels = comma) +
  
  # Theme
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 16, color = "black"),
  
    axis.title.y = element_text(size = 22),
    axis.title.x = element_text(size = 20, margin = margin(t = 10)),
    strip.text.x = element_text(size = 18),
    legend.title = element_text( size = 20, color = "black"),
    legend.text = element_text( size = 16, color = "black"),
    legend.position = "none",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  ) +
  guides(color = "none", fill = "none", shape = "none")+
coord_cartesian(ylim = c(0, 10000)) 

# Print plot
max_bio_LD

2.3 Specific growth rates using all_easylinear

Specific Growth Rates Biotin trial 1 using all_easylinear:

bio1 <- fluor %>%
  filter(Treatment == "-Biotin", Round %in% c(1) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

bio_lin1 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, bio1, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(bio_lin1, log = "y")

summary(bio_lin1)
## $`1:-Biotin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -0.1426 -0.3358  0.3733  0.3348  0.1439 -0.2512 -0.1225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.66156    0.21185  17.284 1.19e-05 ***
## x            0.20749    0.02938   7.063  0.00088 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3109 on 5 degrees of freedom
## Multiple R-squared:  0.9089, Adjusted R-squared:  0.8907 
## F-statistic: 49.88 on 1 and 5 DF,  p-value: 0.0008799
## 
## 
## $`1:-Biotin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -0.2611  0.1266  0.1952  0.1704 -0.1053 -0.1483  0.0224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.77140    0.13104   28.78 9.49e-07 ***
## x            0.21617    0.01817   11.90 7.40e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1923 on 5 degrees of freedom
## Multiple R-squared:  0.9659, Adjusted R-squared:  0.959 
## F-statistic: 141.5 on 1 and 5 DF,  p-value: 7.396e-05
## 
## 
## $`1:-Biotin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5         6         7 
## -0.018655 -0.066394  0.061972  0.130916 -0.079399 -0.037937  0.009497 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.479439   0.142763   17.37 1.16e-05 ***
## x           0.257843   0.007742   33.30 4.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08194 on 5 degrees of freedom
## Multiple R-squared:  0.9955, Adjusted R-squared:  0.9946 
## F-statistic:  1109 on 1 and 5 DF,  p-value: 4.589e-07
coef(bio_lin1)
##                y0    y0_lm     mumax        lag
## 1:-Biotin:1 33.75 38.92188 0.2074863 -0.6871588
## 1:-Biotin:2 33.46 43.44102 0.2161749 -1.2076036
## 1:-Biotin:3 49.69 11.93456 0.2578427  5.5319204
rsquared(bio_lin1)
## 1:-Biotin:1.r2 1:-Biotin:2.r2 1:-Biotin:3.r2 
##      0.9088927      0.9658730      0.9955119

Specific Growth Rates Biotin trial 2 using all_easylinear

bio2 <- fluor %>%
  filter(Treatment == "-Biotin", Round %in% c(2) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

bio_lin2 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, bio2, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(bio_lin2, log = "y")

summary(bio_lin2)
## $`2:-Biotin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -1.29805  0.90387  0.62634  0.53091 -0.40115 -0.02684 -0.33508 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.63417    1.14283   1.430  0.21213   
## x            0.37606    0.07849   4.791  0.00492 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8307 on 5 degrees of freedom
## Multiple R-squared:  0.8211, Adjusted R-squared:  0.7854 
## F-statistic: 22.96 on 1 and 5 DF,  p-value: 0.004922
## 
## 
## $`2:-Biotin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.06137 -0.08479 -0.10821  0.37956 -0.09211  0.27054 -0.30362 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.36851    0.26561  12.682 5.42e-05 ***
## x            0.24466    0.02466   9.921 0.000178 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.261 on 5 degrees of freedom
## Multiple R-squared:  0.9517, Adjusted R-squared:  0.942 
## F-statistic: 98.42 on 1 and 5 DF,  p-value: 0.0001776
## 
## 
## $`2:-Biotin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.24068 -0.16024 -0.25709 -0.01833  0.16293  0.11462 -0.08257 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.10851    0.16756   24.52 2.10e-06 ***
## x            0.29515    0.01873   15.76 1.87e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1983 on 5 degrees of freedom
## Multiple R-squared:  0.9803, Adjusted R-squared:  0.9763 
## F-statistic: 248.2 on 1 and 5 DF,  p-value: 1.873e-05
coef(bio_lin2)
##                 y0     y0_lm     mumax      lag
## 2:-Biotin:1  98.46  5.125203 0.3760631 7.859002
## 2:-Biotin:2  96.62 29.035114 0.2446627 4.914031
## 2:-Biotin:3 107.23 60.856127 0.2951547 1.919209
rsquared(bio_lin2)
## 2:-Biotin:1.r2 2:-Biotin:2.r2 2:-Biotin:3.r2 
##      0.8211458      0.9516543      0.9802552

Specific Growth Rates Biotin trial 3 using all_easylinear:

bio3 <- fluor %>%
  filter(Treatment == "-Biotin", Round %in% c(3) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

bio_lin3 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, bio3, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(bio_lin3, log = "y")

summary(bio_lin3)
## $`3:-Biotin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.07856 -0.22592  0.15259  0.26659  0.21258 -0.23432 -0.09296 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.81040    0.35557   7.904 0.000522 ***
## x            0.23963    0.02156  11.115 0.000103 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2282 on 5 degrees of freedom
## Multiple R-squared:  0.9611, Adjusted R-squared:  0.9533 
## F-statistic: 123.5 on 1 and 5 DF,  p-value: 0.0001028
## 
## 
## $`3:-Biotin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##          1          2          3          4          5          6          7 
##  0.2353339 -0.1112728 -0.0556536 -0.2798767  0.0004855  0.2056337  0.0053499 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.96414    0.23391   12.67 5.44e-05 ***
## x            0.26925    0.01849   14.56 2.76e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1957 on 5 degrees of freedom
## Multiple R-squared:  0.977,  Adjusted R-squared:  0.9723 
## F-statistic:   212 on 1 and 5 DF,  p-value: 2.76e-05
## 
## 
## $`3:-Biotin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.01417 -0.22670  0.78580  0.49481 -1.90075  0.22234  0.61033 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.56645    2.28724   0.248    0.814
## x            0.18591    0.09401   1.978    0.105
## 
## Residual standard error: 0.9949 on 5 degrees of freedom
## Multiple R-squared:  0.4389, Adjusted R-squared:  0.3267 
## F-statistic: 3.911 on 1 and 5 DF,  p-value: 0.1049
coef(bio_lin3)
##                y0     y0_lm     mumax       lag
## 3:-Biotin:1 98.82 16.616595 0.2396343  7.440080
## 3:-Biotin:2 63.40 19.378045 0.2692450  4.402396
## 3:-Biotin:3 82.43  1.761997 0.1859051 20.685290
rsquared(bio_lin3)
## 3:-Biotin:1.r2 3:-Biotin:2.r2 3:-Biotin:3.r2 
##      0.9611016      0.9769567      0.4388915

Specific Growth Rates Biotin trial 4 using all_easylinear

bio4 <- fluor %>%
  filter(Treatment == "-Biotin", Round %in% c(4) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

bio_lin4 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, bio4, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(bio_lin4, log = "y")

summary(bio_lin4)
## $`4:-Biotin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.30806  0.04318 -0.14397 -0.50617 -0.10602  0.24216  0.16276 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.26650    0.53423   2.371 0.063898 .  
## x            0.23271    0.02897   8.032 0.000484 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3066 on 5 degrees of freedom
## Multiple R-squared:  0.9281, Adjusted R-squared:  0.9137 
## F-statistic: 64.51 on 1 and 5 DF,  p-value: 0.0004838
## 
## 
## $`4:-Biotin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.19446  0.54358 -0.06320 -0.20050 -0.37665  0.05642  0.23480 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.80368    0.59222   3.046  0.02857 * 
## x            0.21685    0.03212   6.752  0.00108 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3399 on 5 degrees of freedom
## Multiple R-squared:  0.9012, Adjusted R-squared:  0.8814 
## F-statistic: 45.58 on 1 and 5 DF,  p-value: 0.001082
## 
## 
## $`4:-Biotin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.02364  0.05129  0.03152 -0.14500 -0.02402 -0.04135  0.10391 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.488502   0.120401   20.67 4.91e-06 ***
## x           0.233579   0.008269   28.25 1.04e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08751 on 5 degrees of freedom
## Multiple R-squared:  0.9938, Adjusted R-squared:  0.9925 
## F-statistic: 797.9 on 1 and 5 DF,  p-value: 1.041e-06
coef(bio_lin4)
##                y0     y0_lm     mumax       lag
## 4:-Biotin:1 57.52  3.548409 0.2327066 11.970583
## 4:-Biotin:2 67.67  6.071954 0.2168469 11.118270
## 4:-Biotin:3 97.55 12.043218 0.2335787  8.955711
rsquared(bio_lin4)
## 4:-Biotin:1.r2 4:-Biotin:2.r2 4:-Biotin:3.r2 
##      0.9280693      0.9011555      0.9937725

Specific Growth Rates Thiamine trial 1 using all_easylinear:

thia1 <- fluor %>%
  filter(Treatment == "-Thiamine", Round %in% c(1) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

thia_lin1 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, thia1, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(thia_lin1, log = "y")

summary(thia_lin1)
## $`1:-Thiamine:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.30173 -0.51432 -0.04376  0.21123  0.09790  0.10677 -0.15955 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4499     0.4668   7.391 0.000713 ***
## x             0.2417     0.0283   8.539 0.000363 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2995 on 5 degrees of freedom
## Multiple R-squared:  0.9358, Adjusted R-squared:  0.923 
## F-statistic: 72.91 on 1 and 5 DF,  p-value: 0.0003627
## 
## 
## $`1:-Thiamine:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.31890  0.23129  0.13999 -0.02326  0.25846 -0.25014 -0.03743 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.14430    0.25181   12.49 5.84e-05 ***
## x            0.26660    0.02338   11.40 9.08e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2474 on 5 degrees of freedom
## Multiple R-squared:  0.963,  Adjusted R-squared:  0.9556 
## F-statistic:   130 on 1 and 5 DF,  p-value: 9.079e-05
## 
## 
## $`1:-Thiamine:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.25586  0.36063  0.14809 -0.29385  0.02147 -0.02179  0.04130 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.33255    0.25195   13.23 4.41e-05 ***
## x            0.26006    0.02339   11.12 0.000103 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2476 on 5 degrees of freedom
## Multiple R-squared:  0.9611, Adjusted R-squared:  0.9533 
## F-statistic: 123.6 on 1 and 5 DF,  p-value: 0.0001027
coef(thia_lin1)
##                   y0    y0_lm     mumax       lag
## 1:-Thiamine:1 173.20 31.49697 0.2416543 7.0536952
## 1:-Thiamine:2  44.15 23.20350 0.2666038 2.4129054
## 1:-Thiamine:3  33.19 28.00966 0.2600576 0.6525443
rsquared(thia_lin1)
## 1:-Thiamine:1.r2 1:-Thiamine:2.r2 1:-Thiamine:3.r2 
##        0.9358222        0.9629722        0.9611167

Specific Growth Rates Thiamine trial 2 using all_easylinear

thia2 <- fluor %>%
  filter(Treatment == "-Thiamine", Round %in% c(2) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

thia_lin2 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, thia2, h = 9, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(thia_lin2, log = "y")

summary(thia_lin2)
## $`2:-Thiamine:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4806 -0.2537 -0.1203  0.2472  0.5625 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.63297    0.26430   17.53 4.84e-07 ***
## x            0.24127    0.02348   10.27 1.79e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3638 on 7 degrees of freedom
## Multiple R-squared:  0.9378, Adjusted R-squared:  0.9289 
## F-statistic: 105.5 on 1 and 7 DF,  p-value: 1.79e-05
## 
## 
## $`2:-Thiamine:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35944 -0.06839  0.00212  0.08203  0.24602 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.80494    0.11101   43.28 9.17e-10 ***
## x            0.22326    0.01166   19.15 2.64e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1806 on 7 degrees of freedom
## Multiple R-squared:  0.9813, Adjusted R-squared:  0.9786 
## F-statistic: 366.7 on 1 and 7 DF,  p-value: 2.636e-07
## 
## 
## $`2:-Thiamine:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42503 -0.12665 -0.00432  0.16254  0.40225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.54004    0.16251   27.94 1.93e-08 ***
## x            0.25614    0.01707   15.01 1.40e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2644 on 7 degrees of freedom
## Multiple R-squared:  0.9699, Adjusted R-squared:  0.9656 
## F-statistic: 225.2 on 1 and 7 DF,  p-value: 1.4e-06
coef(thia_lin2)
##                   y0     y0_lm     mumax       lag
## 2:-Thiamine:1 156.08 102.81917 0.2412673 1.7300188
## 2:-Thiamine:2 131.83 122.11205 0.2232630 0.3429772
## 2:-Thiamine:3 140.09  93.69448 0.2561405 1.5704110
rsquared(thia_lin2)
## 2:-Thiamine:1.r2 2:-Thiamine:2.r2 2:-Thiamine:3.r2 
##        0.9378050        0.9812701        0.9698585

Specific Growth Rates Thiamine trial 3 using all_easylinear:

thia3 <- fluor %>%
  filter(Treatment == "-Thiamine", Round %in% c(3) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

thia_lin3 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, thia3, h = 9, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(thia_lin3, log = "y")

summary(thia_lin3)
## $`3:-Thiamine:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35661 -0.21370 -0.07224  0.18665  0.48952 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.79058    0.36271   4.937  0.00168 ** 
## x            0.25375    0.01937  13.101 3.52e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3001 on 7 degrees of freedom
## Multiple R-squared:  0.9608, Adjusted R-squared:  0.9552 
## F-statistic: 171.6 on 1 and 7 DF,  p-value: 3.521e-06
## 
## 
## $`3:-Thiamine:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49865  0.01803  0.07384  0.10434  0.40557 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.81256    0.35153   5.156  0.00132 ** 
## x            0.31014    0.02091  14.833 1.52e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3239 on 7 degrees of freedom
## Multiple R-squared:  0.9692, Adjusted R-squared:  0.9648 
## F-statistic:   220 on 1 and 7 DF,  p-value: 1.517e-06
## 
## 
## $`3:-Thiamine:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4806 -0.1946 -0.0558  0.3160  0.5117 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.99093    0.40297   4.941  0.00167 ** 
## x            0.25389    0.02397  10.593 1.46e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3713 on 7 degrees of freedom
## Multiple R-squared:  0.9413, Adjusted R-squared:  0.9329 
## F-statistic: 112.2 on 1 and 7 DF,  p-value: 1.462e-05
coef(thia_lin3)
##                  y0    y0_lm     mumax      lag
## 3:-Thiamine:1 49.92 5.992945 0.2537505 8.354027
## 3:-Thiamine:2 72.28 6.126119 0.3101377 7.957711
## 3:-Thiamine:3 57.00 7.322315 0.2538908 8.082708
rsquared(thia_lin3)
## 3:-Thiamine:1.r2 3:-Thiamine:2.r2 3:-Thiamine:3.r2 
##        0.9608132        0.9691648        0.9412777

Specific Growth Rates Thiamine trial 4 using all_easylinear

thia4 <- fluor %>%
  filter(Treatment == "-Thiamine", Round %in% c(4) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

thia_lin4 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, thia4, h = 9, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(thia_lin4, log = "y")

summary(thia_lin4)
## $`4:-Thiamine:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34954 -0.14112  0.08196  0.17775  0.22950 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4689     0.1816   19.11 2.68e-07 ***
## x             0.2533     0.0139   18.22 3.71e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2153 on 7 degrees of freedom
## Multiple R-squared:  0.9794, Adjusted R-squared:  0.9764 
## F-statistic:   332 on 1 and 7 DF,  p-value: 3.71e-07
## 
## 
## $`4:-Thiamine:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.269152 -0.100510  0.007589  0.106550  0.281775 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.87449    0.21124   13.61 2.72e-06 ***
## x            0.24309    0.01256   19.35 2.46e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1946 on 7 degrees of freedom
## Multiple R-squared:  0.9816, Adjusted R-squared:  0.979 
## F-statistic: 374.3 on 1 and 7 DF,  p-value: 2.457e-07
## 
## 
## $`4:-Thiamine:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7143  0.0431  0.1993  0.3349  0.4359 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  2.65854    0.76560   3.472   0.0104 * 
## x            0.23859    0.04554   5.239   0.0012 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7055 on 7 degrees of freedom
## Multiple R-squared:  0.7968, Adjusted R-squared:  0.7678 
## F-statistic: 27.45 on 1 and 7 DF,  p-value: 0.0012
coef(thia_lin4)
##                   y0    y0_lm     mumax      lag
## 4:-Thiamine:1 161.92 32.10143 0.2532588 6.389518
## 4:-Thiamine:2  76.34 17.71644 0.2430897 6.008910
## 4:-Thiamine:3 109.65 14.27538 0.2385910 8.544988
rsquared(thia_lin4)
## 4:-Thiamine:1.r2 4:-Thiamine:2.r2 4:-Thiamine:3.r2 
##        0.9793538        0.9816432        0.7968191

Specific Growth Rates Replete trial 1 using all_easylinear:

rep1 <- fluor %>%
  filter(Treatment == "Replete", Round %in% c(1)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

rep_lin1 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, rep1, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(rep_lin1, log = "y")

summary(rep_lin1)
## $`1:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.23658  0.13397  0.19451 -0.01732 -0.04291  0.10941 -0.14107 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.29023    0.23491   14.01 3.34e-05 ***
## x            0.21573    0.01613   13.37 4.19e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1707 on 5 degrees of freedom
## Multiple R-squared:  0.9728, Adjusted R-squared:  0.9674 
## F-statistic: 178.8 on 1 and 5 DF,  p-value: 4.186e-05
## 
## 
## $`1:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.01471  0.09138 -0.04521 -0.16653  0.07629  0.15920 -0.10042 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.03875    0.17358   17.51 1.12e-05 ***
## x            0.22585    0.01192   18.95 7.55e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1262 on 5 degrees of freedom
## Multiple R-squared:  0.9863, Adjusted R-squared:  0.9835 
## F-statistic: 358.9 on 1 and 5 DF,  p-value: 7.551e-06
## 
## 
## $`1:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5         6         7 
##  0.025927  0.009719  0.219121 -0.438641  0.040203  0.154873 -0.011203 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.29711    0.31788   10.37 0.000143 ***
## x            0.22019    0.02183   10.09 0.000164 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.231 on 5 degrees of freedom
## Multiple R-squared:  0.9531, Adjusted R-squared:  0.9438 
## F-statistic: 101.7 on 1 and 5 DF,  p-value: 0.0001641
coef(rep_lin1)
##                y0    y0_lm     mumax      lag
## 1:Replete:1 66.01 26.84907 0.2157266 4.169977
## 1:Replete:2 87.82 20.87905 0.2258538 6.360502
## 1:Replete:3 94.03 27.03446 0.2201909 5.661004
rsquared(rep_lin1)
## 1:Replete:1.r2 1:Replete:2.r2 1:Replete:3.r2 
##      0.9727947      0.9862601      0.9531494

Specific Growth Rates Replete trial 2 using all_easylinear:

rep2 <- fluor %>%
  filter(Treatment == "Replete", Round %in% c(2)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

rep_lin2 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, rep2, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(rep_lin2, log = "y")

summary(rep_lin2)
## $`2:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.18529 -0.03780 -0.16478 -0.18114  0.17936 -0.07894  0.09801 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.37352    0.11553   37.86 2.42e-07 ***
## x            0.18726    0.01602   11.69 8.05e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1695 on 5 degrees of freedom
## Multiple R-squared:  0.9647, Adjusted R-squared:  0.9576 
## F-statistic: 136.6 on 1 and 5 DF,  p-value: 8.053e-05
## 
## 
## $`2:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.04068  0.56895 -0.48634 -0.31279 -0.05303  0.38914 -0.06524 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.36733    0.27564  15.844 1.82e-05 ***
## x            0.20689    0.03822   5.412  0.00291 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4045 on 5 degrees of freedom
## Multiple R-squared:  0.8542, Adjusted R-squared:  0.825 
## F-statistic: 29.29 on 1 and 5 DF,  p-value: 0.002913
## 
## 
## $`2:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.17429  0.11485  0.12945  0.04594 -0.14406  0.10395 -0.07585 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.38836    0.22011   15.39  2.1e-05 ***
## x            0.24033    0.01335   18.01  9.7e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1412 on 5 degrees of freedom
## Multiple R-squared:  0.9848, Adjusted R-squared:  0.9818 
## F-statistic: 324.3 on 1 and 5 DF,  p-value: 9.701e-06
coef(rep_lin2)
##                y0    y0_lm     mumax        lag
## 2:Replete:1 95.47 79.32272 0.1872621  0.9894552
## 2:Replete:2 75.69 78.83269 0.2068881 -0.1966365
## 2:Replete:3 51.77 29.61731 0.2403281  2.3237062
rsquared(rep_lin2)
## 2:Replete:1.r2 2:Replete:2.r2 2:Replete:3.r2 
##      0.9646965      0.8542034      0.9848147

Specific Growth Rates Replete trial 3 using all_easylinear

rep3 <- fluor %>%
  filter(Treatment == "Replete", Round %in% c(3)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

rep_lin3 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, rep3, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(rep_lin3, log = "y")

summary(rep_lin3)
## $`3:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5         6         7 
##  0.056275  0.286364 -0.388386  0.001474 -0.003777 -0.212798  0.260848 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.26798    0.31594   7.179 0.000816 ***
## x            0.26516    0.02498  10.616 0.000128 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2643 on 5 degrees of freedom
## Multiple R-squared:  0.9575, Adjusted R-squared:  0.949 
## F-statistic: 112.7 on 1 and 5 DF,  p-value: 0.0001282
## 
## 
## $`3:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.01109 -0.07478 -0.04567  0.22757 -0.09197 -0.00871 -0.01753 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.69027    0.20383   8.292 0.000416 ***
## x            0.25682    0.01105  23.232 2.75e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.117 on 5 degrees of freedom
## Multiple R-squared:  0.9908, Adjusted R-squared:  0.989 
## F-statistic: 539.7 on 1 and 5 DF,  p-value: 2.75e-06
## 
## 
## $`3:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.37438 -0.38671 -0.12656 -0.27670  0.47488  0.07385 -0.13314 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.09308    0.42459    4.93 0.004361 ** 
## x            0.30144    0.03357    8.98 0.000286 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3552 on 5 degrees of freedom
## Multiple R-squared:  0.9416, Adjusted R-squared:  0.9299 
## F-statistic: 80.65 on 1 and 5 DF,  p-value: 0.0002856
coef(rep_lin3)
##                y0    y0_lm     mumax       lag
## 3:Replete:1 46.52 9.659839 0.2651609  5.928118
## 3:Replete:2 77.27 5.420959 0.2568156 10.346075
## 3:Replete:3 61.72 8.109824 0.3014418  6.732749
rsquared(rep_lin3)
## 3:Replete:1.r2 3:Replete:2.r2 3:Replete:3.r2 
##      0.9575206      0.9908211      0.9416198

Specific Growth Rates Replete trial 4 using all_easylinear:

rep4 <- fluor %>%
  filter(Treatment == "Replete", Round %in% c(4)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

rep_lin4 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, rep4, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(rep_lin4, log = "y")

summary(rep_lin4)
## $`4:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.09195  0.19036  0.07953 -0.20968 -0.15320  0.21726 -0.03232 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.00947    0.25121   11.98 7.15e-05 ***
## x            0.21058    0.01725   12.21 6.53e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1826 on 5 degrees of freedom
## Multiple R-squared:  0.9675, Adjusted R-squared:  0.961 
## F-statistic:   149 on 1 and 5 DF,  p-value: 6.529e-05
## 
## 
## $`4:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.11124  0.04917 -0.34079  0.14249 -0.10230  0.22701 -0.08682 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.08172    0.28910   10.66 0.000126 ***
## x            0.20366    0.01986   10.26 0.000151 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2101 on 5 degrees of freedom
## Multiple R-squared:  0.9546, Adjusted R-squared:  0.9456 
## F-statistic: 105.2 on 1 and 5 DF,  p-value: 0.0001513
## 
## 
## $`4:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.14336  0.08611 -0.35495  0.21552 -0.30373  0.09004  0.12366 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.07289    0.43805   4.732 0.005187 ** 
## x            0.17495    0.02376   7.364 0.000725 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2514 on 5 degrees of freedom
## Multiple R-squared:  0.9156, Adjusted R-squared:  0.8987 
## F-statistic: 54.23 on 1 and 5 DF,  p-value: 0.0007253
coef(rep_lin4)
##                 y0     y0_lm     mumax       lag
## 4:Replete:1  78.16 20.276716 0.2105809  6.407442
## 4:Replete:2 157.37 21.795906 0.2036571  9.706892
## 4:Replete:3  87.34  7.947773 0.1749477 13.700762
rsquared(rep_lin4)
## 4:Replete:1.r2 4:Replete:2.r2 4:Replete:3.r2 
##      0.9675256      0.9546300      0.9155837

Combine all specific growth rate values from each treatment:

extract_results <- function(model, treatment_name) {
  coefs <- coef(model)
  r2vals <- rsquared(model)
  
  df <- as.data.frame(coefs) %>%
    tibble::rownames_to_column("Group") %>%
    mutate(
      Treatment = treatment_name,
      r2 = r2vals
    ) %>%
    # Split "Group" into Round and Rep (Group looks like "1:-Thiamine:2")
    tidyr::separate(Group, into = c("Round", "Treatment_label", "Rep"), sep = ":", remove = FALSE) %>%
    mutate(Round = as.integer(Round))
  
  return(df)
}

# --- Extract Biotin rounds ---
df_bio1 <- extract_results(bio_lin1, "-Biotin")
df_bio2 <- extract_results(bio_lin2, "-Biotin")
df_bio3 <- extract_results(bio_lin3, "-Biotin")
df_bio4 <- extract_results(bio_lin4, "-Biotin")
df_bio <- bind_rows(df_bio1, df_bio2, df_bio3, df_bio4)

# --- Extract Thiamine rounds ---
df_thia1 <- extract_results(thia_lin1, "-Thiamine")
df_thia2 <- extract_results(thia_lin2, "-Thiamine")
df_thia3 <- extract_results(thia_lin3, "-Thiamine")
df_thia4 <- extract_results(thia_lin4, "-Thiamine")
df_thia <- bind_rows(df_thia1, df_thia2, df_thia3, df_thia4)

# --- Extract Replete rounds ---
df_rep1 <- extract_results(rep_lin1, "Replete")
df_rep2 <- extract_results(rep_lin2, "Replete")
df_rep3 <- extract_results(rep_lin3, "Replete")
df_rep4 <- extract_results(rep_lin4, "Replete")
df_rep <- bind_rows(df_rep1, df_rep2, df_rep3, df_rep4)

# --- Combine all treatments ---
all_lin_results <- bind_rows(df_bio, df_thia, df_rep)

# Optional: ensure Round is treated as a factor
all_lin_results$Round <- factor(all_lin_results$Round)

Run an ANOVA on specific growth rate values and Tukey HSD post-hoc test for treatment:

all_lin_results$Round <- factor(all_lin_results$Round)
linear <- aov(mumax ~ Treatment , data = all_lin_results)
summary(linear)
##             Df  Sum Sq  Mean Sq F value Pr(>F)
## Treatment    2 0.00517 0.002586   1.861  0.172
## Residuals   33 0.04587 0.001390
TukeyHSD(linear)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mumax ~ Treatment, data = all_lin_results)
## 
## $Treatment
##                           diff         lwr         upr     p adj
## -Thiamine--Biotin  0.005533672 -0.03181359 0.042880937 0.9298833
## Replete--Biotin   -0.022203946 -0.05955121 0.015143319 0.3234309
## Replete--Thiamine -0.027737618 -0.06508488 0.009609647 0.1780243

Plot specific growth rates for Pyro grown in -Thiamine, -Biotin, and Replete:

all_lin_results$Treatment <- factor(
  all_lin_results$Treatment,
  levels = c("-Cobalamin",  "-Thiamine", "-Biotin","Replete")
)

mu <- ggplot(all_lin_results, aes(x = Treatment, y = mumax, color = Treatment)) +
  # Raw points
  geom_point(size = 3, shape = 16, position = position_jitter(width = 0.1)) +
  # Mean ± SE
  stat_summary(fun = mean, geom = "point", shape = 1, size = 5, stroke = 1.2) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  scale_color_manual(values = c(
    "-Cobalamin" = "#102E50",
    "Replete"    = "#F5C45E",
    "-Thiamine"  = "deeppink4",
    "-Biotin"    = "#667028"
  )) +
  scale_x_discrete(labels = c(
    "-Cobalamin" = expression("\u2212B"[12]),
    "Replete"    = "Replete",
    "-Thiamine"  = expression("\u2212Thiamine"),
    "-Biotin"    = expression("\u2212Biotin")
  )) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    panel.border = element_rect(fill = NA, color = "black"),
    axis.title = element_text( size = 20),
    axis.text = element_text( size = 20, color = "black"),
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text = element_text( size = 20, color = "black"),
    legend.title = element_blank(),
    legend.text = element_blank(),
    legend.key = element_blank()
  ) +
    labs(
      x = "Media Treatment",
      y = expression(paste("(", mu, ") RFU" * "\u00B7" * "day"^-1, " "))
    ) + guides(color = "none", fill = "none", shape = "none")

mu

2.4 Figure 3

bo <- 
  LD_OTB_all_growth / 
  (max_bio_LD +
   mu ) +
  plot_layout(guides = "collect") &
  theme(
    legend.position = "bottom",
    panel.spacing = unit(.1, "lines"),
    plot.margin = margin(15, 15, 15, 15),
    axis.title.y = element_text(margin = margin(r = 12)),
    axis.title.x = element_text(margin = margin(t = 8)),
    legend.margin = margin(t = 4, b = 4),
    legend.box.margin = margin(5, 5, 5, 5),
    axis.text.x = element_text(size = 16),
    
  )
bo

3 P. bahamense grown in five B12 concentrations

Chlorophyll-a fluorescence (RFU) data from P. bahamense cultures were used to visualize growth patterns, quantify maximum biomass, and estimate specific growth rates under media conditions deficient in media lacking one targeted B vitamin. Specific growth rates were computed using the “alleasy_linear” function of the “growthrates” package in R Studio (Hall et al., 2014; RStudio Team, 2020).

Import and Parse the Dataset:

b12 <-as.data.frame(read_csv("b12_levels.csv")) 
# Remove periods and convert to numeric
b12$Time <- as.numeric(sub("T", "", b12$Time))
# Convert Date column to Date type
b12$Date <- mdy(b12$Date)
b12$Treatment <- unlist(b12$Treatment)
#covert timepoints to hours
b12$Time <- b12$Time * 48    # hours
b12$Time <- b12$Time / 24

# put treatment level in desired order
b12 <- b12 %>%
  mutate(Treatment = factor(Treatment, levels = c(
    "369pM",    # Place the levels in desired order
    "36pM",
    "3.69pM",
    "0.369pM",
    "0.0369pM"
  )))

Prepare data for plotting:

levels_order <- c("0.0369pM", "0.369pM", "3.69pM", "36pM", "369pM")
b12$Treatment <- factor(b12$Treatment, levels = levels_order)

# Named color palette
pp_palette <- c(
  "0.0369pM" = "lightblue",
  "0.369pM"  = "#99c2e6",
  "3.69pM"   = "#66a3d2",
  "36pM"     = "#3380bf",
  "369pM"    = "#004080"
)

3.1 Visulaize growth

growth <- ggplot(b12, aes(x = Time, y = `Chlorophyll-A (RFU)`, color = Treatment)) +
  geom_point() +
  geom_smooth(aes(color = Treatment), se = FALSE) +
  facet_grid(Round ~ Treatment) +
  labs(
    x = "Time (days)",
    y = expression(paste("Chlorophyll-", italic("a"), " fluorescence (RFU)"))
  ) +
  
 
   scale_y_continuous(labels = scales::comma) +


  scale_color_manual(values = pp_palette) +
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(size = 12, color = "black"),
    axis.title.y = element_text(size = 20, margin = margin(t = 10)),
    axis.title.x = element_text(size = 20, margin = margin(t = 10)),
    strip.text.x = element_text(size = 16),
    strip.text = element_text( size = 16, color = "black"),
    legend.title = element_text( size = 16, color = "black"),
    legend.text = element_text( size = 16, color = "black"),
    legend.position = "bottom",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt"),
    panel.spacing = unit(.15, "cm")
  )

growth

3.2 Calculate maximum RFU values for each round and medium treatment

max_bio <- b12  %>%
  filter(Treatment %in% c(Treatment, levels = c("369pM", "36pM", "3.69pM", "0.369pM", "0.0369pM"))) %>%
  group_by(Treatment, Bio_Rep) %>% 
  summarise(Max_Chl = max(`Chlorophyll-A (RFU)`, na.rm = TRUE)) %>%
  ungroup()

print(max_bio)
## # A tibble: 15 × 3
##    Treatment Bio_Rep Max_Chl
##    <fct>       <dbl>   <dbl>
##  1 0.0369pM        1    842.
##  2 0.0369pM        2    541.
##  3 0.0369pM        3    928.
##  4 0.369pM         1   1968.
##  5 0.369pM         2   1897.
##  6 0.369pM         3   1866.
##  7 3.69pM          1   1923.
##  8 3.69pM          2   3929.
##  9 3.69pM          3   2146.
## 10 36pM            1   5580.
## 11 36pM            2   6041.
## 12 36pM            3   4801.
## 13 369pM           1   5287.
## 14 369pM           2   4424.
## 15 369pM           3   5539.

Run an ANOVA followed by a Tukey HSD test for the treatment factor

result<- aov(Max_Chl~Treatment, max_bio)
summary(result)
##             Df   Sum Sq  Mean Sq F value   Pr(>F)    
## Treatment    4 49705730 12426432   31.27 1.25e-05 ***
## Residuals   10  3973244   397324                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(result)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Max_Chl ~ Treatment, data = max_bio)
## 
## $Treatment
##                       diff        lwr      upr     p adj
## 0.369pM-0.0369pM 1139.8033  -554.0107 2833.617 0.2494204
## 3.69pM-0.0369pM  1895.5967   201.7826 3589.411 0.0272120
## 36pM-0.0369pM    4703.7467  3009.9326 6397.561 0.0000277
## 369pM-0.0369pM   4312.9433  2619.1293 6006.757 0.0000597
## 3.69pM-0.369pM    755.7933  -938.0207 2449.607 0.6023509
## 36pM-0.369pM     3563.9433  1870.1293 5257.757 0.0003042
## 369pM-0.369pM    3173.1400  1479.3259 4866.954 0.0007806
## 36pM-3.69pM      2808.1500  1114.3359 4501.964 0.0020036
## 369pM-3.69pM     2417.3467   723.5326 4111.161 0.0058828
## 369pM-36pM       -390.8033 -2084.6174 1303.011 0.9366508

Plot maximum biomass:

# Define factor levels (low to high B12)
levels_order <- c("0.0369pM", "0.369pM", "3.69pM", "36pM", "369pM")
max_bio$Treatment <- factor(trimws(max_bio$Treatment), levels = levels_order)

# Calculate summary statistics
summary_stats <- max_bio %>%
  group_by(Treatment) %>%
  summarise(
    mean_chl = mean(Max_Chl, na.rm = TRUE),
    sd_chl   = sd(Max_Chl, na.rm = TRUE),
    .groups = "drop"
  )
summary_stats$Treatment <- factor(summary_stats$Treatment, levels = levels_order)

# Define color palette
pp_palette <- c(
  "0.0369pM" = "lightblue",
  "0.369pM"  = "#99c2e6",
  "3.69pM"   = "#66a3d2",
  "36pM"     = "#3380bf",
  "369pM"    = "#004080"
)




max_bio_plot <- ggplot(summary_stats, aes(x = Treatment)) +
  
  # Replicate points
  geom_point(
    data = max_bio,
    aes(y = Max_Chl, color = Treatment),
    size = 3,
    alpha = 1,
    position = position_jitter(width = 0.1)
  ) +
  
 stat_summary(
  data = max_bio,
  aes(y = Max_Chl, color = Treatment),
  fun = mean,
  geom = "point",
  shape = 1,
  size = 5,
  stroke = 1.2
) +

stat_summary(
  data = max_bio,
  aes(y = Max_Chl, color = Treatment),
  fun.data = mean_sdl,
  fun.args = list(mult = 1),
  geom = "errorbar",
  width = 0.2
)+
  
  # Keep your palette but remove the scale (so no legend)
  scale_color_manual(values = pp_palette, guide = "none") +
  
  scale_y_continuous(labels = scales::comma) +
  
  coord_cartesian(ylim = c(0, 10000)) +
  
  ylab(expression(paste("Maximum RFU"))) +
  xlab(expression(B[12]~Concentration~"(pM)")) +
  
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(size = 20, color = "black", angle = 45, hjust = 1),
    axis.title.y = element_text(size = 20),
    axis.title.x = element_text(size = 20),
    legend.position = "none",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  )

max_bio_plot

3.3 Specific growth rates using all_easylinear

Remove outliers:

# ️⃣ Compute T0 baseline per replicate, treatment, and round
baseline <- aggregate(
  `Chlorophyll-A (RFU)` ~ Bio_Rep + Treatment + Round,
  data = subset(b12, Time == 0),
  FUN = mean,
  na.rm = TRUE
)
names(baseline)[names(baseline) == "Chlorophyll-A (RFU)"] <- "Chl_T0"

# Merge baseline back into the full dataset
b12_merged <- merge(b12, baseline, by = c("Bio_Rep", "Treatment", "Round"), all.x = TRUE)

#  Keep all T0 points, and remove any post-T0 points where replicate < its T0
b12_clean <- subset(
  b12_merged,
  Time == 0 | `Chlorophyll-A (RFU)` >= Chl_T0
)

Specific growth rate with all easy linear for 0.0369PM:

pM_1 <- b12_clean %>%
  filter(
    Treatment %in% c("0.0369pM"),
    Round %in% c(1)
  ) %>%                     
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

mu1 <- all_easylinear(Chl ~ Time | Round + Treatment + Bio_Rep, pM_1, h =4, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(mu1, log = "y")

summary(mu1)
## $`1:0.0369pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
##  0.21354 -0.37851  0.11642  0.04856 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  4.72824    0.52579   8.993   0.0121 *
## x            0.19590    0.07155   2.738   0.1115  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.32 on 2 degrees of freedom
## Multiple R-squared:  0.7894, Adjusted R-squared:  0.6841 
## F-statistic: 7.496 on 1 and 2 DF,  p-value: 0.1115
## 
## 
## $`1:0.0369pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.10998  0.29391 -0.05595 -0.12799 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  5.39799    0.22595  23.890  0.00175 **
## x            0.07310    0.02382   3.069  0.09177 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2429 on 2 degrees of freedom
## Multiple R-squared:  0.8249, Adjusted R-squared:  0.7373 
## F-statistic: 9.421 on 1 and 2 DF,  p-value: 0.09177
## 
## 
## $`1:0.0369pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
##  0.29094 -0.55228  0.23175  0.02959 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   4.5299     0.5772   7.848   0.0159 *
## x             0.2125     0.1054   2.017   0.1813  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4713 on 2 degrees of freedom
## Multiple R-squared:  0.6703, Adjusted R-squared:  0.5055 
## F-statistic: 4.067 on 1 and 2 DF,  p-value: 0.1813
coef(mu1)
##                  y0     y0_lm      mumax        lag
## 1:0.0369pM:1 145.02 113.09593 0.19590188  1.2691828
## 1:0.0369pM:2 197.95 220.96199 0.07310236 -1.5044141
## 1:0.0369pM:3  97.44  92.75275 0.21252290  0.2319725
rsquared(mu1)
## 1:0.0369pM:1.r2 1:0.0369pM:2.r2 1:0.0369pM:3.r2 
##       0.7893916       0.8248814       0.6703408
#not significant slopes

Specific growth rate with all easy linear for 0.369pM:

pM_2 <- b12_clean %>%
  filter(Treatment %in% c("0.369pM"),
         Round %in% c(1),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

mu2 <-all_easylinear(Chl ~ Time|Treatment+Bio_Rep, pM_2, h = 4, quota = 1)


par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(mu2, log = "y")   # line only


summary(mu2)
## $`0.369pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.03871  0.14944 -0.14406  0.03332 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.72896    0.22995  20.566  0.00236 **
## x            0.23253    0.02555   9.101  0.01186 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1512 on 2 degrees of freedom
## Multiple R-squared:  0.9764, Adjusted R-squared:  0.9646 
## F-statistic: 82.83 on 1 and 2 DF,  p-value: 0.01186
## 
## 
## $`0.369pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.03225 -0.15070  0.39814 -0.21519 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  5.12741    0.70102   7.314   0.0182 *
## x            0.20224    0.07559   2.675   0.1159  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3381 on 2 degrees of freedom
## Multiple R-squared:  0.7816, Adjusted R-squared:  0.6724 
## F-statistic: 7.157 on 1 and 2 DF,  p-value: 0.1159
## 
## 
## $`0.369pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.26163  0.40004 -0.01518 -0.12322 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  4.21953    0.72414   5.827   0.0282 *
## x            0.28628    0.07809   3.666   0.0670 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3492 on 2 degrees of freedom
## Multiple R-squared:  0.8705, Adjusted R-squared:  0.8057 
## F-statistic: 13.44 on 1 and 2 DF,  p-value: 0.06701
coef(mu2)
##               y0     y0_lm     mumax      lag
## 0.369pM:1 248.73 113.17754 0.2325283 3.386298
## 0.369pM:2 375.27 168.57963 0.2022351 3.956966
## 0.369pM:3 182.06  68.00143 0.2862760 3.440063
rsquared(mu2)
## 0.369pM:1.r2 0.369pM:2.r2 0.369pM:3.r2 
##    0.9764234    0.7815961    0.8704739
#not significant slopes

Specific growth rate with all easy linear for 3.69PM:

pM_3 <- b12_clean %>%
  filter(Treatment %in% c("3.69pM"),
         Round %in% c(1),
         Time >= 0) %>%  
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)


mu3 <-all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, pM_3, h =5, quota = 1)


par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(mu3 , log = "y")

summary(mu3)
## $`1:3.69pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.03126 -0.21999  0.32699  0.13101 -0.20676 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.60153    0.36017  12.776  0.00103 **
## x            0.24850    0.04245   5.854  0.00994 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2685 on 3 degrees of freedom
## Multiple R-squared:  0.9195, Adjusted R-squared:  0.8927 
## F-statistic: 34.28 on 1 and 3 DF,  p-value: 0.009935
## 
## 
## $`1:3.69pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5 
##  0.100883 -0.081917 -0.005946 -0.145891  0.132871 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.40262    0.10570   51.11 1.65e-05 ***
## x            0.21586    0.02158   10.01  0.00213 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1365 on 3 degrees of freedom
## Multiple R-squared:  0.9709, Adjusted R-squared:  0.9612 
## F-statistic: 100.1 on 1 and 3 DF,  p-value: 0.002126
## 
## 
## $`1:3.69pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5 
## -0.076428  0.082299  0.033975 -0.009133 -0.030712 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.8058     0.1154   41.65 3.05e-05 ***
## x             0.2069     0.0111   18.63 0.000337 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07023 on 3 degrees of freedom
## Multiple R-squared:  0.9914, Adjusted R-squared:  0.9886 
## F-statistic: 347.1 on 1 and 3 DF,  p-value: 0.0003375
coef(mu3)
##                y0     y0_lm     mumax       lag
## 1:3.69pM:1 142.16  99.63649 0.2485048 1.4302532
## 1:3.69pM:2 245.55 221.98664 0.2158566 0.4673628
## 1:3.69pM:3 271.24 122.21605 0.2068770 3.8535635
rsquared(mu3)
## 1:3.69pM:1.r2 1:3.69pM:2.r2 1:3.69pM:3.r2 
##     0.9195173     0.9708995     0.9914315

Specific growth rate with all easy linear for 36PM:

pM_4 <- b12_clean %>% 
  filter(Treatment %in% c("36pM"),
         Round %in% c(1),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)


mu4 <-all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, pM_4, h = 5, quota = 1)


par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(mu4 , log = "y")

summary(mu4)
## $`1:36pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.02088 -0.04810  0.12932 -0.19786  0.09576 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.69402    0.24700   19.00 0.000318 ***
## x            0.23918    0.02377   10.06 0.002090 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1503 on 3 degrees of freedom
## Multiple R-squared:  0.9712, Adjusted R-squared:  0.9616 
## F-statistic: 101.3 on 1 and 3 DF,  p-value: 0.00209
## 
## 
## $`1:36pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5 
## -0.5827  0.4830  0.2809  0.1390 -0.3201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.92981    0.81583   6.043  0.00909 **
## x            0.18334    0.06122   2.995  0.05791 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5078 on 3 degrees of freedom
## Multiple R-squared:  0.7494, Adjusted R-squared:  0.6658 
## F-statistic: 8.969 on 1 and 3 DF,  p-value: 0.05791
## 
## 
## $`1:36pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.05819  0.08475 -0.32792  0.16882  0.01616 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.22326    0.36356  11.617  0.00137 **
## x            0.27122    0.03498   7.753  0.00446 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2213 on 3 degrees of freedom
## Multiple R-squared:  0.9525, Adjusted R-squared:  0.9366 
## F-statistic: 60.11 on 1 and 3 DF,  p-value: 0.004463
coef(mu4)
##              y0     y0_lm     mumax      lag
## 1:36pM:1 330.46 109.29127 0.2391784 4.626125
## 1:36pM:2 205.58 138.35322 0.1833369 2.160096
## 1:36pM:3 104.26  68.25567 0.2712184 1.561941
rsquared(mu4)
## 1:36pM:1.r2 1:36pM:2.r2 1:36pM:3.r2 
##   0.9712271   0.7493526   0.9524610

Specific growth rate with all easy linear for 369PM:

pM_5 <- b12_clean %>% 
  filter(Treatment %in% c("369pM"),
         Round %in% c(1),
         Time >= 0) %>% 
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

mu5 <-all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, pM_5, h = 5, quota = 1)


par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(mu5 , log = "y")

summary(mu5)
## $`1:369pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.08633  0.17016 -0.04861  0.01842 -0.05363 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.74188    0.15680   30.24 7.94e-05 ***
## x            0.25385    0.01538   16.51 0.000484 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1183 on 3 degrees of freedom
## Multiple R-squared:  0.9891, Adjusted R-squared:  0.9855 
## F-statistic: 272.6 on 1 and 3 DF,  p-value: 0.0004836
## 
## 
## $`1:369pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5 
##  0.221300 -0.358207  0.001927 -0.035733  0.170713 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.6394     0.3488  13.303 0.000918 ***
## x             0.2400     0.0342   7.017 0.005946 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2631 on 3 degrees of freedom
## Multiple R-squared:  0.9426, Adjusted R-squared:  0.9234 
## F-statistic: 49.23 on 1 and 3 DF,  p-value: 0.005946
## 
## 
## $`1:369pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.26764  0.51741 -0.14472 -0.19223  0.08718 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.27559    0.49285   8.675  0.00322 **
## x            0.25738    0.05808   4.431  0.02135 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3673 on 3 degrees of freedom
## Multiple R-squared:  0.8675, Adjusted R-squared:  0.8233 
## F-statistic: 19.64 on 1 and 3 DF,  p-value: 0.02135
coef(mu5)
##               y0     y0_lm     mumax      lag
## 1:369pM:1 286.00 114.64947 0.2538508 3.600983
## 1:369pM:2 321.18 103.48174 0.2399582 4.720016
## 1:369pM:3 106.96  71.92225 0.2573820 1.541947
rsquared(mu5)
## 1:369pM:1.r2 1:369pM:2.r2 1:369pM:3.r2 
##    0.9891142    0.9425652    0.8674686

Make dataframe with fits:

  make_df <- function(model, round_num, experiment_name) {
  coefs <- coef(model)
  
  if (is.matrix(coefs)) {
    df <- as.data.frame(coefs) %>%
      tibble::rownames_to_column("Group") %>%
      mutate(
        Round = round_num,
        Experiment = experiment_name,
        r2 = rsquared(model)
      )
  } else {
    df <- data.frame(
      Group = names(coefs),
      Estimate = as.numeric(coefs),
      Round = round_num,
      Experiment = experiment_name,
      r2 = rep(rsquared(model), length(coefs))
    )
  }
  
  return(df)
}

# Example usage for your single fits
df1 <- make_df(mu1, 1, "0.0369pM")
df2 <- make_df(mu2, 1, "0.369pM")
df3 <- make_df(mu3, 1, "3.69pM")
df4 <- make_df(mu4, 1, "36pM")
df5 <- make_df(mu5, 1, "369pM")

# Combine all results
all_results <- bind_rows(df1, df2, df3, df4, df5)

# Check
all_results
##           Group     y0     y0_lm      mumax        lag Round Experiment
## 1  1:0.0369pM:1 145.02 113.09593 0.19590188  1.2691828     1   0.0369pM
## 2  1:0.0369pM:2 197.95 220.96199 0.07310236 -1.5044141     1   0.0369pM
## 3  1:0.0369pM:3  97.44  92.75275 0.21252290  0.2319725     1   0.0369pM
## 4     0.369pM:1 248.73 113.17754 0.23252831  3.3862982     1    0.369pM
## 5     0.369pM:2 375.27 168.57963 0.20223514  3.9569659     1    0.369pM
## 6     0.369pM:3 182.06  68.00143 0.28627601  3.4400634     1    0.369pM
## 7    1:3.69pM:1 142.16  99.63649 0.24850476  1.4302532     1     3.69pM
## 8    1:3.69pM:2 245.55 221.98664 0.21585664  0.4673628     1     3.69pM
## 9    1:3.69pM:3 271.24 122.21605 0.20687699  3.8535635     1     3.69pM
## 10     1:36pM:1 330.46 109.29127 0.23917836  4.6261254     1       36pM
## 11     1:36pM:2 205.58 138.35322 0.18333686  2.1600963     1       36pM
## 12     1:36pM:3 104.26  68.25567 0.27121840  1.5619415     1       36pM
## 13    1:369pM:1 286.00 114.64947 0.25385084  3.6009826     1      369pM
## 14    1:369pM:2 321.18 103.48174 0.23995820  4.7200161     1      369pM
## 15    1:369pM:3 106.96  71.92225 0.25738195  1.5419469     1      369pM
##           r2
## 1  0.7893916
## 2  0.8248814
## 3  0.6703408
## 4  0.9764234
## 5  0.7815961
## 6  0.8704739
## 7  0.9195173
## 8  0.9708995
## 9  0.9914315
## 10 0.9712271
## 11 0.7493526
## 12 0.9524610
## 13 0.9891142
## 14 0.9425652
## 15 0.8674686

ANOVA and post-hoc multiple comparisons of mumax among varying media treatments:

aov_result <- aov(mumax ~ Experiment, data = all_results)
summary(aov_result)
##             Df  Sum Sq  Mean Sq F value Pr(>F)
## Experiment   4 0.01503 0.003757    1.85  0.196
## Residuals   10 0.02031 0.002031
# Post-hoc Tukey test
TukeyHSD(aov_result)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mumax ~ Experiment, data = all_results)
## 
## $Experiment
##                          diff         lwr       upr     p adj
## 0.369pM-0.0369pM  0.079837439 -0.04125361 0.2009285 0.2651098
## 3.69pM-0.0369pM   0.063237082 -0.05785396 0.1843281 0.4652490
## 369pM-0.0369pM    0.089887949 -0.03120310 0.2109790 0.1807721
## 36pM-0.0369pM     0.070735492 -0.05035555 0.1918265 0.3656496
## 3.69pM-0.369pM   -0.016600357 -0.13769140 0.1044907 0.9901101
## 369pM-0.369pM     0.010050511 -0.11104054 0.1311416 0.9985578
## 36pM-0.369pM     -0.009101947 -0.13019299 0.1119891 0.9990215
## 369pM-3.69pM      0.026650867 -0.09444018 0.1477419 0.9459297
## 36pM-3.69pM       0.007498410 -0.11359264 0.1285895 0.9995435
## 36pM-369pM       -0.019152457 -0.14024350 0.1019386 0.9831961

T-test between 0.0369pM and all other concentrations:

all_results$Group <- ifelse(all_results$Experiment == "0.0369pM", 
                            "A", 
                            "B")
all_results$Group <- factor(all_results$Group)

t.test(mumax~ Group, data = all_results)
## 
##  Welch Two Sample t-test
## 
## data:  mumax by Group
## t = -1.6943, df = 2.1577, p-value = 0.2231
## alternative hypothesis: true difference in means between group A and group B is not equal to 0
## 95 percent confidence interval:
##  -0.2558391  0.1039901
## sample estimates:
## mean in group A mean in group B 
##       0.1605090       0.2364335

Plot mumax with all easy linear data:

# Define the order of treatments from least to greatest
levels_order <- c("0.0369pM", "0.369pM", "3.69pM", "36pM", "369pM")

# Convert Treatment column to a factor with this order
all_results$Experiment <- factor(all_results$Experiment, levels = levels_order)

# calculate standard deviation and mean 
b12_summary <- all_results %>%
  group_by(Experiment) %>%
  summarise(
    mean_mumax = mean(mumax, na.rm = TRUE),
    se_mumax   = sd(mumax, na.rm = TRUE)/sqrt(n()),
    .groups = "drop"
  )

mu_b12 <- ggplot(all_results, aes(x = Experiment, y = mumax, color = Experiment)) +
  
  # Raw replicate points (jittered)
  geom_point(size = 3, shape = 16, 
             position = position_jitter(width = 0.1), alpha = 1) +
  
  # Mean per Experiment+Round (open circle)
  stat_summary(fun = mean, geom = "point", shape = 1, size = 5, stroke = 1.2)  +
  
  # Error bars (standard error)
  stat_summary(fun.data = mean_se, geom = "errorbar",
               width = 0.2) +
  
  # Color palette — WITH legend removed
  scale_color_manual(
    values = c(
      "0.0369pM" = "lightblue",
      "0.369pM"  = "#99c2e6",
      "3.69pM"   = "#66a3d2",
      "36pM"     = "#3380bf",
      "369pM"    = "#004080"
    ),
    guide = "none"   # ← THIS removes the legend safely
  ) +

  # Safe visible range
  coord_cartesian(ylim = c(0, 0.4)) +
  
  # Labels + theme
  labs(
    x = expression(B[12]~Concentration~"(pM)"),
    y = expression(paste(mu, " (RFU" %*% "day"^-1, ")"))
  ) +
  
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    panel.border = element_rect(fill = NA, color = "black"),
    axis.title = element_text(size = 20),
    axis.text = element_text(size = 20, color = "black"),
    axis.text.x = element_text(size = 20, color = "black", angle = 45, hjust = 1),
    strip.text = element_text(size = 20, color = "black"),

    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  )

mu_b12

3.4 Figure 4

# Ensure BOTH subplots have no legend inside patchwork
max_bio_nl <- max_bio_plot + theme(legend.position = "none")
mu_nl      <- mu_b12 + theme(legend.position = "none")

# Build the layout
bo <- 
  (
    growth /                
    (max_bio_nl | mu_nl)    
  ) +
  plot_layout(
    widths = c(4, 2),        # your original width ratio
    guides = "collect"        # ← THIS is what actually gathers one legend
  ) &
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    panel.spacing = unit(0.4, "lines"),
    plot.margin = margin(15, 15, 15, 15),
    axis.title.y = element_text(margin = margin(r = 12)),
    axis.title.x = element_text(margin = margin(t = 8)),
    legend.margin = margin(t = 4, b = 4),
    legend.box.margin = margin(5, 5, 5, 5),
    legend.ticks = element_blank(),
    axis.text.x = element_text(size = 16)
  )

bo

4 P. bahamense recruitment experiments

Low-diveristy P. bahamense culture was co-cultured with bacteria recruited from OTB (Old Tampa Bay) water, Indian River Lagoon (IRL) culture, and IRL surface seawater and grown in replete media and media missing B12

Import data:

OTB water:

OTB <- as.data.frame(read_csv("OTB_1006.csv"))

OTB$Time <- as.numeric(sub('.', '', OTB$Time))
OTB$Date <- ymd(OTB$Date)
OTB$Time <- OTB$Time * 48

OTB %<>%
  filter(Treatment %in% c("-Cobalamin", "Replete")) %>%
  mutate(Treatment = factor(Treatment, levels = c("-Cobalamin", "Replete")))

OTB$Time_in_Days <- OTB$Time / 24

IRL water:

IRLrec<-as.data.frame(read_csv("IRL_1006.csv"))
IRLrec$Time <- as.numeric(sub('.', '', IRLrec$Time))
IRLrec$Date <- ymd(IRLrec$Date)
IRLrec$Time <- IRLrec$Time * 48
IRLrec$Time_in_Days <- IRLrec$Time / 24  # Convert Time from hours to days

IRL culture:

IRLcul <- as.data.frame(read_csv("1003_1006.csv"))
IRLcul$Time <- as.numeric(gsub("\\.", "", IRLcul$Time))  # Remove all periods
IRLcul$Date <- ymd(IRLcul$Date)
IRLcul$Time <- IRLcul$Time * 48
IRLcul$Time_in_Days <- IRLcul$Time / 24 

Combine all datasets into one with a new ‘Dataset’ column:

all_experiments <- bind_rows(
  fluor %>% mutate(Dataset = "Low Diversity OTB"),
  IRLcul %>% mutate(Dataset = "IRL culture microbiome"),
  OTB %>% mutate(Dataset = "OTB water microbiome"),
  IRLrec %>% mutate(Dataset = "IRL water microbiome")
) %>%mutate(
  Treatment = recode(Treatment,
                     "-Cobalamin" = "\u2013B12",
                     "-Biotin" = "\u2013Biotin",
                     "-Thiamine" = "\u2013Thiamine"),
  Treatment = factor(Treatment, levels = c("\u2013B12", "Replete"))
) %>%
  filter(Round %in% c(1, 2), Treatment %in% c("\u2013B12", "Replete"))

4.1 Visulaize growth of P. bahamense with recruited bacteria

# Filter out "Low Diversity OTB"
all_experiments_filtered <- all_experiments %>%
  filter(Dataset != "Low Diversity OTB") %>%
  mutate(
    Treatment = recode(Treatment,
                       "–B12" = "—B12",  # convert en dash to em dash
                       "Replete" = "Replete"),
    Treatment = factor(Treatment, levels = c("—B12", "Replete")),
    Dataset = factor(Dataset, levels = c(
      "IRL culture microbiome",
      "OTB water microbiome",
      "IRL water microbiome"
    ))
  )

growth <- ggplot(all_experiments_filtered, aes(x = Time_in_Days, y = `Chlorophyll-A (RFU)`, color = Treatment)) +
  geom_point(size = 3) +
  geom_smooth(se = FALSE) +
  scale_color_manual(values = c("—B12" = "#102E50", "Replete" = "#F5C45E")) +
  facet_grid(Round ~ Dataset) +
  labs(
    x = "Time (days)",
    y = expression(paste("Chlorophyll-", italic("a"), " fluorescence (RFU)"))
  ) +
  scale_y_continuous(labels = label_comma()) +
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 16, color = "black"),
    axis.text.x = element_text(size = 12, color = "black"),
    axis.title.y = element_text(size = 20, margin = margin(t = 10)),
    axis.title.x = element_text(size = 20, margin = margin(t = 10)),
    strip.text.x = element_text(size = 20),
    strip.text = element_text(size = 20, color = "black"),
    legend.title = element_text(size = 16, color = "black"),
    legend.text = element_text(size = 16, color = "black"),
    legend.position = "bottom",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt"),
    panel.spacing = unit(.1, "cm") ) +
  coord_cartesian(ylim = c(0, 5000))

plot(growth)

Create new variable for experimemt and combine experimental data into dataframe:

fluor$Experiment <- "fluor"  
IRLcul$Experiment <- "IRLcul"  
OTB$Experiment <- "OTB"  
IRLrec$Experiment <- "IRLwat"  
all_exp <- bind_rows(fluor, IRLcul, OTB, IRLrec) %>%filter(Round %in% c(1, 2, 3,4)) %>%
  filter(Treatment %in% c("-Cobalamin", "Replete")) 

4.2 Calculate maximum RFU values for each round and medium treatment

max_bio <- all_exp %>%
  group_by(Treatment,  Bio_Rep,Experiment, Round) %>%
  summarise(Max_Chl = max(`Chlorophyll-A (RFU)`, na.rm = TRUE), .groups = "drop") %>%
  arrange(Treatment,  Bio_Rep)
print(max_bio)
## # A tibble: 44 × 5
##    Treatment  Bio_Rep Experiment Round Max_Chl
##    <chr>        <dbl> <chr>      <dbl>   <dbl>
##  1 -Cobalamin       1 IRLcul         1   103. 
##  2 -Cobalamin       1 IRLcul         2   136. 
##  3 -Cobalamin       1 IRLwat         1   893. 
##  4 -Cobalamin       1 IRLwat         2   904. 
##  5 -Cobalamin       1 OTB            1   341. 
##  6 -Cobalamin       1 OTB            2    87.8
##  7 -Cobalamin       2 IRLcul         1   147. 
##  8 -Cobalamin       2 IRLcul         2   173. 
##  9 -Cobalamin       2 IRLwat         1   674. 
## 10 -Cobalamin       2 IRLwat         2  1126. 
## # ℹ 34 more rows

Plot maximum biomass

# Define experiment order
experiment_order <- c("fluor", "IRLcul", "OTB", "IRLwat")

# Filter and prepare data for -Cobalamin only
plot_data <- max_bio %>%
  filter(Treatment == "-Cobalamin") %>%
  mutate(Experiment = factor(Experiment, levels = experiment_order))
# 
# plot_summary <- summary_stats %>%
#   filter(Treatment == "-Cobalamin") %>%
#   mutate(Experiment = factor(Experiment, levels = experiment_order))

# Build plot
max_bio_plot <- ggplot(plot_data, aes(x = Experiment, y = Max_Chl, color = Treatment)) +

  # Raw replicate points (solid)
  geom_point(size = 5, shape = 16, alpha = 1,
             position = position_jitter(width = 0.15)) +

  stat_summary(
    fun = mean,
    geom = "point",
    shape = 1,
    size = 5,
    stroke = 1.2
  ) +
  
  # Error bars: mean ± SD  (NOT SE)
  stat_summary(
    fun.data = mean_sdl,
    fun.args = list(mult = 1),   # 1 SD above/below mean
    geom = "errorbar",
    width = 0.2
  ) +

  # Legend with open circle
  scale_color_manual(
    values = c("-Cobalamin" = "#102E50"),
    labels = c("-Cobalamin" = expression("—B"[12])),
    name = "Treatment",
    guide = guide_legend(
      override.aes = list(
        shape = 21,
        fill = "white",   # open circle in legend
        stroke = 1.2,
        size = 5
      )
    )
  ) +

  # Axis labels
  labs(
    x = "Source",
    y = expression(paste("Maximum RFU"))
  ) +

  # X-axis labels
  scale_x_discrete(labels = c(
    fluor = "Low Diversity",
    IRLcul = "IRL Culture",
    OTB = "OTB Seawater",
    IRLwat = "IRL Seawater"
  )) +

  # Y-axis formatting
  scale_y_continuous(labels = label_comma()) +

  # # Facet by Round
  # facet_grid(~Round) +

  # Theme
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 14, color = "black"),
    axis.title.y = element_text(size = 20),
    axis.title.x = element_text(size = 20, margin = margin(t = 10)),
    strip.text.x = element_text(size = 20),
    strip.text = element_text(size = 20, color = "black"),
    legend.title = element_text(size = 16, color = "black"),
    legend.text = element_text(size = 16, color = "black"),
    legend.position = "bottom",
    legend.key.size = unit(1.5, "cm"),
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  ) +

  # Set visible y-axis range safely
  coord_cartesian(ylim = c(50, 1500))

# Print plot
max_bio_plot

Run an ANOVA followed by a Tukey HSD test for the experiment factor:

max_B12 <- all_exp %>%
  group_by(Experiment ,Treatment, Round, Bio_Rep) %>% 
  summarise(Max_Chl = max(`Chlorophyll-A (RFU)`, na.rm = TRUE)) %>%
  ungroup()  %>% filter(Treatment %in% c("-Cobalamin"))


aov_biomass <- aov(Max_Chl ~ Experiment, data =max_B12)
summary(aov_biomass)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## Experiment   3 2481763  827254   38.95 4.45e-08 ***
## Residuals   18  382327   21240                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov_biomass)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Max_Chl ~ Experiment, data = max_B12)
## 
## $Experiment
##                     diff       lwr       upr     p adj
## IRLcul-fluor  -148.55167 -414.4357  117.3324 0.4145013
## IRLwat-fluor   679.19000  413.3060  945.0740 0.0000057
## OTB-fluor      -23.03833 -288.9224  242.8457 0.9946502
## IRLwat-IRLcul  827.74167  589.9278 1065.5556 0.0000001
## OTB-IRLcul     125.51333 -112.3006  363.3272 0.4626003
## OTB-IRLwat    -702.22833 -940.0422 -464.4144 0.0000007

4.3 Specific growth rates using all_easylinear

Remove outliers from IRL culture data: outliers are considered values the drop below T0 RFU values

# Compute T0 baseline per replicate, treatment, and round
baseline <- aggregate(
  `Chlorophyll-A (RFU)` ~ Bio_Rep + Treatment + Round,
  data = subset(IRLcul, Time == 0),
  FUN = mean,
  na.rm = TRUE
)
names(baseline)[names(baseline) == "Chlorophyll-A (RFU)"] <- "Chl_T0"

# Merge baseline back into the full dataset
IRLcul_merged <- merge(IRLcul, baseline, by = c("Bio_Rep", "Treatment", "Round"), all.x = TRUE)

# Keep all T0 points, and remove any post-T0 points where replicate < its T0
IRLcul_clean <- subset(
  IRLcul_merged,
  Time == 0 | `Chlorophyll-A (RFU)` >= Chl_T0
)

IRL culture trial 1 Replete:

IRLcul_mumax <- IRLcul_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(1),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)
IRLcul1 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRLcul_mumax, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLcul1 , log = "y")


summary(IRLcul1 )
## $`1:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.03861  0.16034 -0.15112  0.16012 -0.18606 -0.07377  0.12911 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.7592193  0.3438351   5.116  0.00372 ** 
## x           0.0075785  0.0006407  11.828  7.6e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1627 on 5 degrees of freedom
## Multiple R-squared:  0.9655, Adjusted R-squared:  0.9586 
## F-statistic: 139.9 on 1 and 5 DF,  p-value: 7.603e-05
## 
## 
## $`1:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.05977  0.09664  0.02457  0.12515 -0.33697  0.07566  0.07473 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.793494   0.469100   1.692 0.151524    
## x           0.007143   0.000691  10.336 0.000146 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1755 on 5 degrees of freedom
## Multiple R-squared:  0.9553, Adjusted R-squared:  0.9463 
## F-statistic: 106.8 on 1 and 5 DF,  p-value: 0.0001459
## 
## 
## $`1:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.03760 -0.17761  0.53725 -0.39227  0.01679  0.10783 -0.05440 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 1.488243   0.895579   1.662  0.15745   
## x           0.005008   0.001233   4.062  0.00971 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3132 on 5 degrees of freedom
## Multiple R-squared:  0.7674, Adjusted R-squared:  0.7209 
## F-statistic:  16.5 on 1 and 5 DF,  p-value: 0.009713
coef(IRLcul1 )
##                y0    y0_lm       mumax      lag
## 1:Replete:1 74.22 5.807901 0.007578469 336.1912
## 1:Replete:2 88.61 2.211109 0.007142610 516.7230
## 1:Replete:3 71.52 4.429309 0.005007900 555.4690
rsquared(IRLcul1)
## 1:Replete:1.r2 1:Replete:2.r2 1:Replete:3.r2 
##      0.9654963      0.9552895      0.7674161

IRL culture trial 2 Replete:

IRLcul_mumax2 <- IRLcul_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(2),
         Time >= 0) %>% 
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

IRLcul2 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRLcul_mumax2, h =7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLcul2 , log = "y")


summary(IRLcul2 )
## $`2:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.23499  0.44031 -0.16656 -0.01381 -0.09558  0.10722 -0.03659 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.7712453  0.4712913   5.880 0.002020 ** 
## x           0.0069299  0.0009628   7.198 0.000806 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2445 on 5 degrees of freedom
## Multiple R-squared:  0.912,  Adjusted R-squared:  0.8944 
## F-statistic: 51.81 on 1 and 5 DF,  p-value: 0.0008063
## 
## 
## $`2:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.15234 -0.55703  0.04284  0.55888  0.12964 -0.23617 -0.09050 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 1.626180   0.806476   2.016  0.09983 . 
## x           0.008378   0.001503   5.575  0.00256 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3817 on 5 degrees of freedom
## Multiple R-squared:  0.8614, Adjusted R-squared:  0.8337 
## F-statistic: 31.08 on 1 and 5 DF,  p-value: 0.002558
## 
## 
## $`2:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.31177 -0.51835 -0.13281  0.31474  0.18935 -0.07054 -0.09416 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 2.906698   0.628814   4.623  0.00572 **
## x           0.006738   0.001285   5.245  0.00334 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3263 on 5 degrees of freedom
## Multiple R-squared:  0.8462, Adjusted R-squared:  0.8155 
## F-statistic: 27.51 on 1 and 5 DF,  p-value: 0.003339
coef(IRLcul2 )
##                y0     y0_lm       mumax      lag
## 2:Replete:1 73.64 15.978520 0.006929859 220.4869
## 2:Replete:2 81.69  5.084417 0.008378133 331.4284
## 2:Replete:3 71.42 18.296276 0.006738242 202.1121
rsquared(IRLcul2)
## 2:Replete:1.r2 2:Replete:2.r2 2:Replete:3.r2 
##      0.9119822      0.8614255      0.8462232

Remove outliers from OTB water data:

# Compute T0 baseline per replicate, treatment, and round
baseline <- aggregate(
  `Chlorophyll-A (RFU)` ~ Bio_Rep + Treatment + Round,
  data = subset(OTB, Time == 0),
  FUN = mean,
  na.rm = TRUE
)
names(baseline)[names(baseline) == "Chlorophyll-A (RFU)"] <- "Chl_T0"

# Merge baseline back into the full dataset
OTBwat_merged <- merge(OTB, baseline, by = c("Bio_Rep", "Treatment", "Round"), all.x = TRUE)

# Keep all T0 points, and remove any post-T0 points where replicate < its T0
OTBwat_clean <- subset(
  OTBwat_merged,
  Time == 0 | `Chlorophyll-A (RFU)` >= Chl_T0
)

OTB water trial 1 Replete:

OTB_mumax1 <- OTBwat_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(1),
         Time >= 0) %>% 
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)


OTBwat1 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, OTB_mumax1, h = 7, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(OTBwat1 , log = "y")


summary(OTBwat1)
## $`1:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.15859  0.22093 -0.08966 -0.03358  0.16311 -0.01998 -0.08224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.332309   0.210366   20.59 5.00e-06 ***
## x           0.007191   0.000602   11.95 7.25e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1529 on 5 degrees of freedom
## Multiple R-squared:  0.9661, Adjusted R-squared:  0.9594 
## F-statistic: 142.7 on 1 and 5 DF,  p-value: 7.251e-05
## 
## 
## $`1:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.33247 -0.31312 -0.07991 -0.10915  0.06674  0.08438  0.01858 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.6685656  0.3405662  10.772 0.000120 ***
## x           0.0084842  0.0008604   9.861 0.000183 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2185 on 5 degrees of freedom
## Multiple R-squared:  0.9511, Adjusted R-squared:  0.9413 
## F-statistic: 97.23 on 1 and 5 DF,  p-value: 0.0001828
## 
## 
## $`1:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.05109  0.05233 -0.25963  0.01366  0.20820  0.01297 -0.07862 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.3310777  0.2440626   17.75 1.04e-05 ***
## x           0.0068159  0.0006166   11.05 0.000106 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1566 on 5 degrees of freedom
## Multiple R-squared:  0.9607, Adjusted R-squared:  0.9528 
## F-statistic: 122.2 on 1 and 5 DF,  p-value: 0.0001055
coef(OTBwat1)
##                 y0    y0_lm       mumax      lag
## 1:Replete:1 211.13 76.11986 0.007190613 141.8745
## 1:Replete:2 101.61 39.19564 0.008484247 112.2759
## 1:Replete:3 328.76 76.02618 0.006815913 214.8282
rsquared(OTBwat1)
## 1:Replete:1.r2 1:Replete:2.r2 1:Replete:3.r2 
##      0.9661413      0.9510923      0.9606888

OTB water trial 2 Replete:

OTB_mumax2 <- OTBwat_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(2),
         Time >= 0) %>% 
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)


OTBwat2 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, OTB_mumax2, h = 5, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(OTBwat2  , log = "y")



summary(OTBwat2)
## $`2:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.04205 -0.14349  0.12612  0.01002 -0.03470 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.0096618  0.0889572  56.315 1.23e-05 ***
## x           0.0041960  0.0007566   5.546   0.0116 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1148 on 3 degrees of freedom
## Multiple R-squared:  0.9111, Adjusted R-squared:  0.8815 
## F-statistic: 30.76 on 1 and 3 DF,  p-value: 0.01156
## 
## 
## $`2:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.12058  0.06702 -0.35136  0.01932  0.14443 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.724234   0.180952  26.108 0.000123 ***
## x           0.007660   0.001539   4.977 0.015586 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2336 on 3 degrees of freedom
## Multiple R-squared:  0.892,  Adjusted R-squared:  0.856 
## F-statistic: 24.77 on 1 and 3 DF,  p-value: 0.01559
## 
## 
## $`2:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5 
##  0.1065  0.1344 -0.4463  0.0634  0.1420 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.409273   0.224859  19.609  0.00029 ***
## x           0.007442   0.001912   3.891  0.03009 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2903 on 3 degrees of freedom
## Multiple R-squared:  0.8347, Adjusted R-squared:  0.7795 
## F-statistic: 15.14 on 1 and 3 DF,  p-value: 0.03009
coef(OTBwat2)
##                 y0     y0_lm       mumax      lag
## 2:Replete:1 156.29 149.85405 0.004195997 10.02180
## 2:Replete:2 127.08 112.64419 0.007660107 15.74165
## 2:Replete:3  91.45  82.20965 0.007442250 14.31284
rsquared(OTBwat2)
## 2:Replete:1.r2 2:Replete:2.r2 2:Replete:3.r2 
##      0.9111292      0.8919806      0.8346501

Remove outliers from IRL water data:

#Compute T0 baseline per replicate, treatment, and round
baseline <- aggregate(
  `Chlorophyll-A (RFU)` ~ Bio_Rep + Treatment + Round,
  data = subset(IRLrec, Time == 0),
  FUN = mean,
  na.rm = TRUE
)
names(baseline)[names(baseline) == "Chlorophyll-A (RFU)"] <- "Chl_T0"

# Merge baseline back into the full dataset
IRLwat_merged <- merge(IRLrec, baseline, by = c("Bio_Rep", "Treatment", "Round"), all.x = TRUE)

# Keep all T0 points, and remove any post-T0 points where replicate < its T0
IRLwat_clean <- subset(
  IRLwat_merged,
  Time == 0 | `Chlorophyll-A (RFU)` >= Chl_T0
)

IRL water trial 1 Replete:

IRL_mumax1 <- IRLwat_clean %>%
  filter(Treatment %in% c("Replete"), Round %in% c(1)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

IRLwat1 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRL_mumax1, h = 5, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLwat1 , log = "y")




summary(IRLwat1)
## $`1:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.07587  0.09142 -0.13804 -0.22946  0.35195 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.005693   0.268748  18.626 0.000338 ***
## x           0.008304   0.001006   8.258 0.003718 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2644 on 3 degrees of freedom
## Multiple R-squared:  0.9579, Adjusted R-squared:  0.9438 
## F-statistic:  68.2 on 1 and 3 DF,  p-value: 0.003718
## 
## 
## $`1:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.01569 -0.18763  0.22688 -0.29711  0.27355 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.6061505  0.2231989   25.12 0.000138 ***
## x           0.0064232  0.0008246    7.79 0.004403 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2887 on 3 degrees of freedom
## Multiple R-squared:  0.9529, Adjusted R-squared:  0.9372 
## F-statistic: 60.68 on 1 and 3 DF,  p-value: 0.004403
## 
## 
## $`1:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.32328  0.13706  0.35205 -0.08458 -0.08125 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.230102   0.228469  18.515 0.000344 ***
## x           0.013307   0.001943   6.848 0.006374 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.295 on 3 degrees of freedom
## Multiple R-squared:  0.9399, Adjusted R-squared:  0.9198 
## F-statistic: 46.89 on 1 and 3 DF,  p-value: 0.006374
coef(IRLwat1)
##                 y0     y0_lm       mumax        lag
## 1:Replete:1 108.14 149.26049 0.008304417 -38.806618
## 1:Replete:2 267.86 272.09478 0.006423229  -2.442075
## 1:Replete:3 178.44  68.72424 0.013306580  71.705157
rsquared(IRLwat1)
## 1:Replete:1.r2 1:Replete:2.r2 1:Replete:3.r2 
##      0.9578639      0.9528881      0.9398718

IRL water trial 2 Replete:

IRL_mumax2 <- IRLwat_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(2),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

IRLwat2 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRL_mumax2, h =6, quota = 1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLwat2 , log = "y")

summary(IRLwat2)
## $`2:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6 
##  0.27656 -0.39465 -0.12399  0.08819  0.39137 -0.23748 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.126211   0.317395  13.000 0.000202 ***
## x           0.009866   0.001698   5.811 0.004365 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3409 on 4 degrees of freedom
## Multiple R-squared:  0.8941, Adjusted R-squared:  0.8676 
## F-statistic: 33.76 on 1 and 4 DF,  p-value: 0.004365
## 
## 
## $`2:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6 
##  0.03290  0.04031 -0.06486  0.15651 -0.44416  0.27930 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.883134   0.200405  24.366 1.68e-05 ***
## x           0.006915   0.001379   5.014  0.00741 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2769 on 4 degrees of freedom
## Multiple R-squared:  0.8628, Adjusted R-squared:  0.8284 
## F-statistic: 25.14 on 1 and 4 DF,  p-value: 0.007414
## 
## 
## $`2:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6 
## -0.10408  0.24742  0.10827 -0.59427  0.39445 -0.05179 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.516376   0.359111  12.577  0.00023 ***
## x           0.008128   0.001921   4.231  0.01336 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3857 on 4 degrees of freedom
## Multiple R-squared:  0.8174, Adjusted R-squared:  0.7717 
## F-statistic:  17.9 on 1 and 4 DF,  p-value: 0.01336
coef(IRLwat2)
##                 y0     y0_lm       mumax       lag
## 2:Replete:1 102.18  61.94276 0.009865993 50.732380
## 2:Replete:2 136.46 132.04390 0.006914858  4.757449
## 2:Replete:3 105.76  91.50335 0.008128438 17.813607
rsquared(IRLwat2)
## 2:Replete:1.r2 2:Replete:2.r2 2:Replete:3.r2 
##      0.8940791      0.8627528      0.8173779

IRL water trial 1 -B12:

IRL_mumax1_b <- IRLwat_clean %>%
  filter(Treatment %in% c("-Cobalamin"),
         Round %in% c(1),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

IRLwat1_b <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRL_mumax1_b, h = 4, quota =1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLwat1_b , log = "y")

summary(IRLwat1_b)
## $`1:-Cobalamin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.13839  0.02356  0.36806 -0.25322 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 4.629046   0.277043  16.709  0.00356 **
## x           0.011942   0.003085   3.871  0.06073 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3311 on 2 degrees of freedom
## Multiple R-squared:  0.8822, Adjusted R-squared:  0.8233 
## F-statistic: 14.98 on 1 and 2 DF,  p-value: 0.06073
## 
## 
## $`1:-Cobalamin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
##  0.09071 -0.12902 -0.06650  0.10480 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.1181197  0.1239394  41.295 0.000586 ***
## x           0.0053753  0.0007698   6.982 0.019900 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1419 on 2 degrees of freedom
## Multiple R-squared:  0.9606, Adjusted R-squared:  0.9409 
## F-statistic: 48.76 on 1 and 2 DF,  p-value: 0.0199
## 
## 
## $`1:-Cobalamin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
##  0.06389 -0.33062  0.46957 -0.20284 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 3.332327   0.711542   4.683   0.0427 *
## x           0.015666   0.004035   3.883   0.0604 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.433 on 2 degrees of freedom
## Multiple R-squared:  0.8829, Adjusted R-squared:  0.8243 
## F-statistic: 15.08 on 1 and 2 DF,  p-value: 0.06038
coef(IRLwat1_b)
##                    y0     y0_lm       mumax       lag
## 1:-Cobalamin:1  89.18 102.41631 0.011941636 -11.58880
## 1:-Cobalamin:2 182.88 167.02102 0.005375287  16.87548
## 1:-Cobalamin:3 112.12  28.00342 0.015665946  88.55150
rsquared(IRLwat1_b)
## 1:-Cobalamin:1.r2 1:-Cobalamin:2.r2 1:-Cobalamin:3.r2 
##         0.8822319         0.9605951         0.8828865

IRL water trial 2 -B12:

IRL_mumax2_b <-  IRLwat_clean %>%
  filter(Treatment %in% c("-Cobalamin"),
         Round %in% c(2),
         Time >= 0) %>%  
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

IRLwat2_b <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRL_mumax2_b, h =5, quota =1)

par(mfrow = c(3, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLwat2_b , log = "y")

summary(IRLwat2_b)
## $`2:-Cobalamin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.08915 -0.35337  0.17809  0.34733 -0.26120 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 4.371518   0.460476   9.493  0.00248 **
## x           0.008698   0.002261   3.847  0.03101 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3432 on 3 degrees of freedom
## Multiple R-squared:  0.8314, Adjusted R-squared:  0.7752 
## F-statistic:  14.8 on 1 and 3 DF,  p-value: 0.03101
## 
## 
## $`2:-Cobalamin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.04172  0.31220 -0.17814 -0.37176  0.27941 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.022854   0.290615  17.284 0.000422 ***
## x           0.007184   0.001842   3.899 0.029934 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3402 on 3 degrees of freedom
## Multiple R-squared:  0.8352, Adjusted R-squared:  0.7803 
## F-statistic:  15.2 on 1 and 3 DF,  p-value: 0.02993
## 
## 
## $`2:-Cobalamin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.23747  0.05598 -0.51117 -0.09547  0.31319 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 4.206086   0.396161  10.617  0.00179 **
## x           0.010498   0.002488   4.219  0.02434 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3777 on 3 degrees of freedom
## Multiple R-squared:  0.8558, Adjusted R-squared:  0.8077 
## F-statistic:  17.8 on 1 and 3 DF,  p-value: 0.02434
coef(IRLwat2_b)
##                    y0     y0_lm       mumax       lag
## 2:-Cobalamin:1 107.43  79.16373 0.008697853 35.103053
## 2:-Cobalamin:2 145.64 151.84401 0.007183847 -5.806909
## 2:-Cobalamin:3 130.09  67.09340 0.010498480 63.070164
rsquared(IRLwat2_b)
## 2:-Cobalamin:1.r2 2:-Cobalamin:2.r2 2:-Cobalamin:3.r2 
##         0.8314296         0.8352092         0.8557594

Extract and combine all specific growth rates:

# Helper function to handle matrix or vector coefficients
make_df <- function(model, round_num, experiment_name) {
  coefs <- coef(model)
  
  if (is.matrix(coefs)) {
    df <- as.data.frame(coefs) %>%
      tibble::rownames_to_column("Group") %>%
      mutate(
        Round = round_num,
        Experiment = experiment_name,
        r2 = rsquared(model)
      )
  } else {
    df <- data.frame(
      Group = names(coefs),
      Estimate = as.numeric(coefs),
      Round = round_num,
      Experiment = experiment_name,
      r2 = rep(rsquared(model), length(coefs))
    )
  }
  
  return(df)
}

# Create data frames for all models with experiment names
df1 <- make_df(IRLcul1, 1, "IRL Culture")
df2 <- make_df(IRLcul2, 2, "IRL Culture")
df3 <- make_df(OTBwat1, 1, "OTB Water")
df4 <- make_df(OTBwat2, 2, "OTB Water")
df5 <- make_df(IRLwat1, 1, "IRL Water")
df6 <- make_df(IRLwat2, 2, "IRL Water")
df7 <- make_df(IRLwat1_b, 1, "IRL Water B")
df8 <- make_df(IRLwat2_b, 2, "IRL Water B")

# Combine all results
all_results <- bind_rows(df1, df2, df3, df4, df5, df6, df7, df8)

# Optional: reorder columns
all_results <- all_results %>%
  select(Experiment, Round, Group, y0, mumax, r2)

Prepare IRL seawater data for T-test for both trials in -B12 media:

df5 <- as.data.frame(coef(IRLwat1))
df5$Round <- 1
df5$Experiment <- "IRL Water"

df6 <- as.data.frame(coef(IRLwat2))
df6$Round <- 2
df6$Experiment <- "IRL Water"

df7 <- as.data.frame(coef(IRLwat1_b))
df7$Round <- 1
df7$Experiment <- "IRL Water -B12"

df8 <- as.data.frame(coef(IRLwat2_b))
df8$Round <- 2
df8$Experiment <- "IRL Water -B12"

# Combine
IRLwater <- bind_rows(df5, df6, df7, df8)

T-test between specific growthrates both trials for IRL water in -B12:

t.test(mumax ~ Experiment, data = IRLwater, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  mumax by Experiment
## t = -0.59074, df = 10, p-value = 0.5678
## alternative hypothesis: true difference in means between group IRL Water and group IRL Water -B12 is not equal to 0
## 95 percent confidence interval:
##  -0.005105451  0.002965605
## sample estimates:
##      mean in group IRL Water mean in group IRL Water -B12 
##                  0.008823919                  0.009893842

Plot IRL seawater specific growth rates:

# Extract the columns needed
mumax_results <- all_results %>%
  select(Experiment, Round, Group, mumax) %>%
  mutate(ExperimentRound = paste(Experiment, Round, sep = "_"))

mumax_results <- mumax_results %>%
  mutate(
    Treatment = factor(
      Experiment,
      levels = c("IRL Water B", "IRL Water")  # order: B12 first, Replete second
    )
  )

# Filter IRL Water only
mumax_IRLwater <- mumax_results %>%
  filter(grepl("IRL Water", Experiment))


IRL_mu <- ggplot(mumax_IRLwater, aes(x = Treatment, y = mumax, color = Treatment)) +
  # Raw replicate points
  geom_point(size = 5, shape = 16, position = position_jitter(width = 0.1)) +
  
  # Mean per Experiment+Round (open circle)
  stat_summary(fun = mean, geom = "point", shape = 1, size = 5, stroke = 1.2) +
  
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  
  # facet_grid(~Round) +
  scale_color_manual(values = c(
    "IRL Water B" = "#102E50",
    "IRL Water" = "#F5C45E",
    "-Thiamine" = "deeppink4",
    "-Biotin" = "#667028"
  )) +
  scale_x_discrete(labels = c(
  "IRL Water B"  = expression("\u2013B"[12]), 
  "IRL Water"   = "Replete",
  "-Thiamine"   = expression("\u2013Thiamine"),
  "-Biotin"     = expression("\u2013Biotin")    
))+
  theme_minimal()  +

  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    panel.border = element_rect(fill = NA, color = "black"),
    axis.title = element_text(size = 20),
    axis.text = element_text(size = 20, color = "black"),
    axis.text.x = element_text( size = 18, color = "black", angle = 45, hjust = 1),
    strip.text = element_text( size = 20, color = "black"),
    legend.title = element_blank(),
    legend.text = element_blank(),
    legend.key = element_blank(),
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  ) +
  labs(
    x = "Media Treatment",
    y = expression(paste("(", mu, ") RFU" * "\u00B7" * "day"^-1, " "))
  ) +
  guides(color = "none", fill = "none", shape = "none") +
  
  # Set visible range safely
  coord_cartesian(ylim = c(0, 0.5)) 
  
IRL_mu

4.4 Figure 5

# Ensure BOTH subplots have no legend inside patchwork
max_bio__2 <- max_bio_plot + theme(legend.position = "none")
mu_2      <- IRL_mu + theme(legend.position = "none")

bo <- 
  growth / 
  (max_bio_plot  +
   IRL_mu ) + 
  plot_layout(
    guides = "collect",
    widths = c(4, 2)     
  ) &
  theme(
    legend.position = "bottom",
    panel.spacing = unit(.1, "lines"),
    plot.margin = margin(15, 15, 15, 15),
    axis.title.y = element_text(margin = margin(r = 12)),
    axis.title.x = element_text(margin = margin(t = 8)),
    legend.margin = margin(t = 4, b = 4),
    legend.box.margin = margin(5, 5, 5, 5),
    axis.text.x = element_text(size = 16)
  )

bo

5 Characterizing P. bahamense recruited bacterial communities

  1. Full-length 16S rRNA was sequenced using Oxford Nanopore MinIon v10.4 flow cells (Oxford Nanopore Technology 2008)
  2. Sequence data was processed and classified using Kraken v2.1.2 (Lu and Salzberg 2020)

Import data:

Recruit_16s <- read.csv("/Users/lydiaruggles/Desktop/KRAKEN-16s-counts-genus.csv")
# Identify sample columns
sample_cols <- names(Recruit_16s)[2:(which(names(Recruit_16s) == "total") - 1)]

# Pivot to long data
Recruit_16s_long <- Recruit_16s %>%
  pivot_longer(
    cols = all_of(sample_cols),
    names_to = "sample",
    values_to = "count"
  )

# Separate samples into new columns
Recruit_16s_long <- Recruit_16s_long %>%
  separate(sample, into = c("experiment","treatment","filter_size","bio_rep"),
           sep = "_", remove = FALSE)

# One genus in one sample with the total count of that genus in that sample
grouped_data1 <- Recruit_16s_long %>%
  mutate(count = as.numeric(count)) %>%
  group_by(treatment, experiment, filter_size, bio_rep, genus) %>%
  summarise(total_count = sum(count), .groups = "drop"
)

Prepare data as wide matrix

# Pivot wider
wide_data1 <- grouped_data1 %>%
  pivot_wider(names_from = genus, values_from = total_count, values_fill = 0)

# Create rownames and remove metadata columns for distance matrix
rownames(wide_data1) <- paste(wide_data1$treatment, wide_data1$experiment, wide_data1$filter_size, wide_data1$bio_rep, sep = "_")

wide_matrix1 <- wide_data1[, !(names(wide_data1) %in% c("treatment", "experiment", "filter_size", "bio_rep"))]
wide_matrix1 <- as.data.frame(lapply(wide_matrix1, as.numeric))

wide_matrix1[is.na(wide_matrix1)] <- 0


# Bray-Curtis + PCoA
bray_dist1 <- vegdist(wide_matrix1, method = "bray")
pcoa_result1 <- pcoa(bray_dist1)

# Build plot data
pcoa_data1 <- as.data.frame(pcoa_result1$vectors)
pcoa_data1$treatment <- wide_data1$treatment
pcoa_data1$experiment <- wide_data1$experiment
pcoa_data1$filter_size <- as.factor(wide_data1$filter_size)

Permanova:

#ensure wide matrix is numeric
wide_matrix1 <- as.data.frame(lapply(wide_matrix1, as.numeric))
wide_matrix1[is.na(wide_matrix1)] <- 0

#Build the Bray–Curtis distance matrix
bray_dist <- vegdist(wide_matrix1, method = "bray")

#Bring metadata back in (from wide_data1)
metadata <- wide_data1[, c("treatment", "experiment", "filter_size", "bio_rep")]


permanova_result <- adonis2(
  bray_dist ~ experiment,
  data = metadata,
  permutations = 999
)

print(permanova_result)
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = bray_dist ~ experiment, data = metadata, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)    
## Model     4   6.2820 0.67542 16.127  0.001 ***
## Residual 31   3.0189 0.32458                  
## Total    35   9.3009 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.1 Figure 6 PCoA plot

sw_pcoa <- ggplot(
  pcoa_data1,
  aes(
    x = Axis.1,
    y = Axis.2,
    color = experiment,
    shape = treatment,
    alpha = filter_size   # keep as factor
  )
) +
  geom_point(size = 5) +
  labs(
    x = paste0("PCoA1 (", round(pcoa_result1$values$Relative_eig[1] * 100, 2), "%)"),
    y = paste0("PCoA2 (", round(pcoa_result1$values$Relative_eig[2] * 100, 2), "%)"),
    color = "Experiment",
    shape = "Treatment",
    alpha = "Filter size"
  ) +
  stat_ellipse(aes(group = experiment)) +
  scale_color_hue() +

  # Map discrete alpha manually
  scale_alpha_manual(values = c("02" = 0.4, "10" = 1)) +

  theme_minimal() +
  theme(
    legend.key.size = unit(1.5, "cm"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 20),
    panel.border = element_rect(color = 'black', fill = NA),
    axis.title.x = element_text(size = 20, color = "black"),
    axis.title.y = element_text(size = 20, color = "black"),
    strip.text = element_text(size = 20, color = "black"),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(size = 20, color = "black")
  ) +
  geom_hline(yintercept = 0, linetype = "dotted", linewidth = 0.5) +
  geom_vline(xintercept = 0, linetype = "dotted", linewidth = 0.5) +
  theme(panel.grid = element_blank())

plot(sw_pcoa)

Prepare data for relative abundance plot:

#Standardize experiment names
grouped_data <- grouped_data1 %>%
  mutate(
    experiment = case_when(
      grepl("1003.1006", experiment, ignore.case = TRUE) ~ "1003.1006",
      grepl("OTB", experiment, ignore.case = TRUE) ~ "OTB.1006",
      grepl("IRL", experiment, ignore.case = TRUE) ~ "IRL.1006",
      TRUE ~ experiment
    ),
    experiment = factor(experiment,
                        levels = c("1003.1006", "OTB_seawater", "OTB.1006", "IRL_seawater", "IRL.1006"),
                        labels = c("Microbiome from IRL culture",
                                   "OTB Seawater",
                                   "Microbiome from OTB seawater",
                                   "IRL seawater",
                                   "Microbiome from IRL seawater"))
  )

#Pivot to long format
sample_cols <- names(Recruit_16s)[2:(which(names(Recruit_16s) == "total") - 1)]

Recruit_16s_long <- Recruit_16s %>%
  pivot_longer(
    cols = all_of(sample_cols),
    names_to = "sample",
    values_to = "count"
  ) %>%
  separate(sample, into = c("experiment","treatment","filter_size","bio_rep"),
           sep = "_", remove = FALSE) %>%
  mutate(count = as.numeric(count))  # ensure numeric

# ️ Summarize counts per sample & genus
grouped_data <- Recruit_16s_long %>%
  group_by(experiment, treatment, filter_size, bio_rep, order) %>%
  summarise(total_count = sum(count), .groups = "drop") %>%
  group_by(experiment, treatment, filter_size, bio_rep) %>%
  mutate(RA = (total_count / sum(total_count)) * 100) %>%
  ungroup() %>%
  filter(treatment %in% c("noB12", "replete", "seawater"))

# ⃣Collapse low abundance genera (<1%) into "z<1% abund."
grouped_data$order <- as.character(grouped_data$order)
grouped_data$order [grouped_data$RA < 1] <- "z<1% abund."
grouped_data$order <- factor(grouped_data$order )  # back to factor

#  Optional: fix experiment & filter labels for plotting
grouped_data <- grouped_data %>%
  mutate(
    experiment = factor(experiment,
                        levels = c("X1003.1006","OTB","OTB.1006","IRL","IRL.1006"),
                        labels = c("Recruited IRL culture",
                                   "OTB Seawater",
                                   "Recruited OTB seawater",
                                   "IRL Seawater",
                                   "Recruited IRL seawater")),
    Treatment = recode(treatment, "-B12" = "\u2014B12"),
    filter_size = recode(filter_size, "10" = "10.0um", "02" = "0.2um"),
    filter_size = factor(filter_size, levels = c("10.0um", "0.2um")),  # <-- top first
    x_label = paste(Treatment, "Rep", bio_rep, sep = " ")
  )

#  Define color palette (ensure enough colors for all genera)
unique_genera <- levels(grouped_data$order )
# Use your custom palette if available, otherwise RColorBrewer
library(RColorBrewer)
n_colors <- length(unique_genera)
palette_colors <- brewer.pal(min(12, n_colors), "Set3")  # adjust if >12 genera
if(n_colors > 12) {
  palette_colors <- colorRampPalette(palette_colors)(n_colors)
}
names(palette_colors) <- unique_genera

5.2 Figure 7 Relative abundance plot

taxa_bar_plot <- ggplot(grouped_data,
                        aes(x = x_label, y = RA, fill = order)) +
  geom_bar(stat = "identity", position = "stack") +
  facet_grid(filter_size ~ experiment, scales = "free_x") +
  theme_classic() +
  theme(
    legend.position = "right",
    legend.title = element_blank(),
    axis.text.x = element_text(size = 18, angle = 45, hjust = 1, color = "black"),
    axis.text.y = element_text(size = 18, color = "black"),
    axis.title.y = element_text(size = 18),
    strip.text = element_text(size = 12),
    legend.text = element_text(size = 18),
    strip.position = "top"
  ) +
  scale_y_continuous() +        # Corrected
  scale_fill_viridis_d(option = "turbo") +             # Deep distinct colors
  guides(fill = guide_legend(ncol = 2)) +
  ylab("Relative Abundance (%)") +
  labs(x = NULL)


taxa_bar_plot